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We review recent theoretical work on two closely related issues: excitation of an isolated 
(DJQ' quantum condensed matter system driven adiabatically across a continuous quantum phase 

transition or a gapless phase, and apparent relaxation of an excited system after a sudden 
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in the physics of ultracold atoms with potential applications in the adiabatic quantum state 
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1. Introduction 

Experiments with ultracold atoms provide unprecedented opportunities to emulate 
precisely tailored model Hamiltonians at temperatures close to zero [1]. These 
model quantum systems are not only well isolated from their environment, but also 
their parameters can be changed in time in a prescribed way [2-9]. Consequently, 
the emerging field of non-equilibrium quantum systems driven out of their initial 
ground state by a time-dependent Hamiltonian is no longer of purely academic 
interest. Quite to the contrary, it becomes essential for such practical applications 
as adiabatic quantum state preparation for quantum simulation [1, 9-13] and, in 
particular, adiabatic quantum computation [14]. In this article we review recent 
theoretical work on two closely related issues: excitation of a condensed matter 
system driven across a continuous quantum phase transition or a gapless phase, 
and apparent relaxation of an excited system towards thermal equilibrium or a 
non-thermal steady state. 

The former issue had been intensively studied both theoretically and experimen- 
tally in different systems at finite temperature including cosmological phase tran- 
sitions, liquid crystals, superfluid helium, convection cells, superconductors, and 
Bose-Einstein condensation, where the excitation is described by the Kibble-Zurek 
mechanism (KZM). In this review we collected recent evidence that, although in 
many situations many details are different than in the finite temperature KZM, its 
key ingredients - which are crossovers between adiabatic and non-adiabatic stages 
of time evolution - remain the same in the zero temperature quantum limit. Con- 
sequently, the excitation of a system in general scales with a power of a transition 
rate. We also include examples where the KZM paradigm is not so predictive, like 
when driving the system across a gapless bosonic phase, even though the excitation 
is still a power of the driving rate, and an example of a disordered system where 
it is the KZM that predicts a logarithmic dependence on the rate instead of the 
usual power law. Finally, we mention an example where the evolution is essentially 
non-adiabatic, i.e., the excitation density diverges with increasing system size. 

Once a system got excited the latter issue becomes important: does an excited 
state of an isolated quantum system relax to any steady state? A quick "global" 
answer is "no" , because unitary evolution cannot evolve an initial pure state into a 
mixed state. However, this does not exclude some sort of apparent local relaxation, 
where expectation values of local obscrvablcs relax to become the same as in a 
global mixed state. With this distinction in mind, we can further ask what is the 
nature of the global steady state? In this article we briefly review recent work 
on integrable and non-integrable quantum systems. General arguments together 
with exactly solvable examples suggest that non-interacting quadratic bosonic or 
fermionic Hamiltonians appear to relax locally to a non-thermal generalised Gibbs 
ensemble (GGE). The ensemble is strongly constrained by integrals of motion which 
is this case are the numbers of Bogoliubov quasiparticles. In case of non-integrable 
systems the evidence is less conclusive. When an excited state of a many-body 
system after a sudden quench has a narrow energy distribution then, according 
to the conjectured "eigenstate thermalisation hypothesis" (ETH), for simple few- 
body observables the state appears to relax to a microcanonical ensemble. This 
thermalisation mechanism is qualitatively diflFerent than the ergodicity required of 
classical non-integrable systems. 

Even though in practice the two issues cannot be always sharply "disentangled" , 
because the relaxation begins already during the excitation of a system, this review 
is divided into two corresponding major parts in Sections 2 and 3. 
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2. Dynamics of a quantum phase transition 
2.1. Introduction 

A quantum phase transition is a fundamental change in the ground state of a 

quantum system when one of parameters in its Hamiltonian is driven across a 
critical point [15]. A continuous (or second order) quantum phase transition can 
often be characterised by an energy gap between the ground state and the first 
excited state or, more generally, a relevant energy scale which vanishes at the 
critical point in the thermodynamic limit of infinite system size. Consequently, no 
matter how slowly the parameter is driven the quantum state of the system cannot 
follow adiabatically the instantaneous ground state near the critical point. It is 
the aim of this part to quantify the level of this inevitable excitation in different 
systems. 

This research is motivated in part by adiabatic quantum computation [14] , where 
one would like to prepare a system in a simple ground state of an initial Hamiltonian 
Hi, and then drive the system adiabatically to a final Hamiltonian Hf, whose 
ground state is a solution to a non-trivial computational problem. Unfortunately, 
since Hi and Hf are qualitatively different, they are generally separated by a critical 
point, where the driving cannot be adiabatic in the limit of infinite system size. 
Consequently, there is an upper limit to the size of the system (or number of qubits) 
which can cross the critical point adiabatically at a given transition rate and, 
consequently, the minimal necessary time of "adiabatic compTitation" increases 
with the number of qubits. The same problem can arise in the context of quantum 
simulators [12, 16], where the key idea is to emulate a model condensed matter 
Hamiltonian with, say, ultracold atoms in an optical lattice potential [1, 8, 13], ions 
in an ion trap [9, 10] or NMR molecules [11]. In this context, one usually prepares 
the system in the uncorrelated ground state of a simple initial Hamiltonian H^, 
like e.g. a Bose-Einstein condensate, and then drives it into a correlated state by 
a slow change of a parameter in the Hamiltonian [2, 6, 8, 9, 11]. 

A quantum phase transition across an isolated quantum critical point between 
two gapped phases turns out to be well described by a quantum version of the 
Kibble-Zurek mechanism (KZM). The essence of the mechanism is an adiabatic- 
impulse-adiabatic approximation, where evolution of a system driven slowly across 
a phase transition is assumed adiabatic in the gapped phases sufficiently far from 
criticality, and impulse in a close neighbourhood of the gapless critical point, where 
the state of the system does not change in this approximation. Consequently, 
the state after the transition can be argued to have a finite correlation length 

^ ~ Tq'"' , where tq is a characteristic time of the adiabatic transition between 
the two phases and i^, z are the critical exponents. This universal scale of length 
determines density of excitations and other physical quantities. For example, the 
density of excitations or excitation energy above the ground state scales as an 
inverse power of the transition time tq. 

However, the world of quantum critical phenomena is rich: phase diagrams have 
gapless lines and gapless phases, and a system driven across a gapless line or phase 
also gets excited. It turns out that in many of such non-standard situations the 
adiabatic-impulse-adiabatic approximation still applies and is able to predict the 
correct scaling of the density of excitations (or density of excitation energy) with 
the transition time, but there are also systems where this approximation is not 
practical or not predictive enough. These include some bosonic systems, which 
are non-adiabatic in the thermodynamic limit, i.e., their density of excitations 
diverges with increasing system size. Nevertheless, the adiabatic-impulse-adiabatic 
approximation is the leading motif of Section 2. 
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Section 2 is organised as follows. In Section 2.2 we briefly review KZM in the 
historical context of finite temperature classical phase transitions, and then in 
Section 2.3 we generalise KZM to a zero temperature quantum phase transition 
between two gapped phases. Once KZM is reformulated in this most "canonical" 
set-up its key ingredient, which is the adiabatic-impulse-adiabatic approximation, 
can be bend and twisted in many different ways as required by actual application. 
Some general examples are described in Sections 2.4, 2.6, 2.7, 2.8, and 2.5. In 
Section 2.4 we drive a finite system across the critical value of the parameter in its 
Hamiltonian and estimate how slow the driving rate has to be for the transition 
to be adiabatic. In Section 2.6 we drive an infinite system into a gapped phase 
starting exactly at a critical point. In Section 2.7 a quantum phase transition takes 
place not in time but in space, i.e., the parameter in the Hamiltonian driving the 
transition is time-independent but inhomogeneous in space, so that different parts 
of the system are in different phases. In Section 2.8 the phase transition is driven in 
an inhomogeneous way such that some parts of the system cross the critical point 
earlier than the other. The inhomogeneous transition turns out to be a way to 
suppress excitations at the critical point. Another way to minimise excitations is a 
transition with a non-linear time-dependence of the driving parameter considered 
in Section 2.5. 

After this incomplete review of non-standard generalisations, we return to the 
basics and in Sections 2.9,2.10,2.11 investigate relations between KZM and the 
Landau-Zener (LZ) model of level anti-crossing [17]. We begin in Section 2.9 by 
a review of the standard LZ model together with more exotic transitions that 
begin or end at the anti-crossing centre. In Section 2.10 we apply the adiabatic- 
impulse-adiabatic approximation, which is central to KZM, to obtain approximate 
solutions of the LZ model in some limits of parameters. Finally, in Section 2.3 
we consider a general quadratic fermionic Hamiltonian which can be mapped to 
sets of independent LZ transitions. This class of intcgrablc models provides exact 
solutions supporting KZM or its generalisations. In the following Section 2.12 we 
rederive KZM within the time-dependent perturbation theory. This perturbative 
treatment reproduces approximately the exact results in the integrable models, but 
it does not seem limited to the quadratic Hamiltonians. Section 2.12 completes the 
general part of Section 2, and the following Sections describe applications of the 
general ideas to a number of specific integrable and non-integrable models. 

The list of examples is opened in Section 2.13 by one of the cornerstones of 
the theory of quantum phase transitions - the integrable quantum Ising chain in 
transverse magnetic field. It has a quantum phase transition between two gapped 
phases: paramagnetic and ferromagnetic. Since the model can be mapped via the 
Jordan- Wigner transformation [18] to a quadratic fermionic Hamiltonian, virtually 
all the general ideas mentioned above are illustrated in the subsections of Section 
2.13. This is why the quantum Ising chain takes relatively large portion of this 
article. 

Once the review of the Ising chain is completed, in Sections 2.14 and 2.15 we 
review another two quadratic fermionic systems: the XY spin chain and the 2D 
Kitaev model. They provide exactly solvable examples where KZM requires care- 
ful generalisation to obtain correct scaling of the density of excitations with the 
transition rate. The next example in Section 2.16 is the random quantum Ising 
chain, where the density of excitations turns out to be a logarithmic function of 
the transition rate. There is no usual KZ power-law scaling in this model, but the 
logarithmic dependence can still be obtained from the adiabatic-impulse-adiabatic 
approximation like in KZM. This example is followed by topological insulators in 
Section 2.17, where the edge zero modes that must exist in a topological insulator 



August 10, 2010 0:9 Advances in Physics revised2 



6 Dynamics of a Quantum Phase Transition and Relaxation to a Steady State 

have a different scaling of the excitation probabihty than in the bulk of the same in- 
sulator. In Section 2.18 the quantum KZM is extended to the Lipkin-Meshkov-Glick 
model with infinite coordination number, and in Section 2.19 to the non-integrable 
Bose-Hubbard model. 

The last model is important because it is another cornerstone of the theory 
of quantum phase transitions and it was realised experimentally in Refs. [2, 13, 
19]. The ID version of the model has the Bcrczinski-Kosterlitz-Thoulcss transition 
between the gapped Mott insulator and gapless superfluid. Since the model is 
non-integrable, we review different approximate results for transitions between the 
two phases. Another closely related model is considered in Section 2.20, where a 
quasi- ID gas of bosonic atoms is loaded into an optical lattice potential. The gas 
is a gapless (critical) Luttinger liquid for which a weak optical lattice potential 
is a relevant perturbation leading to an effective sine-Gordon model with a finite 
gap in its excitation spectrum. This is an example of a transition starting at a 
critical point and going into a gapped phase introduced in Section 2.6. In Section 
2.21 we review mean-field treatment of the transition from a paramagnetic to a 
ferromagnetic phase in a spin-1 Bose-Einstein condensate. 

Finally, in Section 2.22 we consider adiabatic sweep across a gapless phase in a 
quadratic model of harmonic oscillators, and in Section 2.23 a semiclassical treat- 
ment of a many-body Landau-Zener model. The ID gapless harmonic model turns 
out to be non-adiabatic in a sense that its density of excitations increases with the 
system size. The many-body LZ model is in turn an example, where the density of 
LZ anti-crossings is so large that they cannot be treated as independent transitions, 
so the whole problem has to be and can be treated by semiclassical methods. A 
summary of Section 2 is made in Section 2.24. 



2.2. KZM in a classical phase transition 

Phase transition is a fundamental change in the state of a system when one of its 
parameters passes through the critical point. In a second order phase transition, 
the fundamental change is continuous and the critical point is characterised by 
divergent correlation length and relaxation time. This critical slowing down implies 
that no matter how slowly a system is driven across the transition, its evolution 
cannot be adiabatic close to the critical point. As a result, ordering of the state after 
a transition from a disordered symmetric phase to an ordered broken symmetry 
phase is not perfect: the state is a mosaic of ordered domains whose finite size 
^ depends on the rate of the transition. This scenario was first described in the 
cosmological context by Kibble [20] who appealed to relativistic causality to set 
the size of the domains. The dynamical mechanism relevant for second order phase 
transitions was proposed by Zurek [21] and it is briefly as follows. 

A dimensionless distance from the critical point at a finite temperature Tc can 
be defined as 



T-T, 



^ = ■ (1) 



Tc 



In the initial symmetric phase above Tc there are finite thermal fluctuations of the 
order parameter with diverging correlation length 



(2) 
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and relaxation time 

r ~ r (3) 

describing reaction time of the system to external perturbations. This divergent 
relaxation time becomes infinitely slow at the critical point. 

The system can be driven below cither directly by cooling, or by increasing Tc 
above a constant T by varying another parameter such as pressure [22, 23]. In any 
case, there is a time-dependent e(t) running continuously from an initial > to 
a final e/ < 0. Near the critical point it can be linearised as 

e{t) «^ - :f , (4) 

where the coefficient tq can be identified as a "quench time" . This linearised form 
implies a relative transition rate = which diverges near the critical 

point. The evolution of the system with a time-dependent e(t) can be divided into 
three stages, see Fig. 1. Initially the transition rate is slower than the relaxation 
time T and the system adiabatically follows the instantaneous state of thermal 
equilibrium for current e{t). This adiabatic stage lasts until 



'Q 



(5) 



when the transition rate equals the instantaneous relaxation rate r^^ of the 
system and long wavelength fluctuations begin to go out of equilibrium. After e 
reactions of the system are too sluggish to follow the varying e{t) and, in a first 
impulse approximation, the state does not change until — e when the reactions of 
the system become faster than the transition rate again. In this way the system 
arrives at — e, which is already in the symmetry broken phase, still remaining in 
the state of thermal equilibrium at the +e in the symmetric phase, where there are 
small thermal fluctuations of the order parameter with a finite correlation length 

I ~ e-'' ~ . (6) 

These small fiuctuations are the initial state for the last adiabatic stage of the 
evolution after — e. 

Around — e the system "realises" at last that zero order parameter is no longer 
its state of equilibrium. The proper equilibrium state has a finite order parameter 
whose eventual variations in space can be characterised by a finite healing length 
I — el*^^, which is incidentally the same as the correlation length ^ (but see Ref. [24]). 
Near — e the small fluctuations of the order parameter with wave lengths longer 
than the healing length ^ blow up quasi-exponentially in time until magnitude of 
the order parameter becomes comparable to its flnite equilibrium value. At this 
point the order parameter becomes a mosaic of ordered domains whose average 
size is set by The orientation of the order parameter is approximately constant 
in each domain, but uncorrelated between different domains. 

In the finite temperature context the focus was on topological defects. When a 
vacuum manifold of the order parameter has a non-trivial homotopy group, then 
there are stable topological defects. Their density after a transition is set by the size 
^ of the correlated domains [20]. For example, the density of point-like monopoles 
in 3D is ~ and vortices in 2D is 2± Topologically stable defects are a 
robust and relatively easy to detect imprint of the non-equilibrium phase transition. 
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— Reaction time x 

— Inverse transition rate le/(de/dt)l 




Distance from the critical point 8 



Figure 1. The reaction time r and the inverse transition rate |e/e| as a function of dimensionless distance 
e from the critical point. The reaction time equals the inverse transition rate at +e and —e. These two 
points mark crossovers between the impulse and adiabatic stages of time evolution. 



However, most of them are not permanent because, for instance, pairs of vortices 
with opposite winding number attract each other and annihilate. This equilibration 
process, known as phase ordering kinetics [25], leads to gradual coarse-graining of 
the order parameter on increasing scale of length. What eventually remains is a 
net winding number of the initial mosaic. This is why in the finite temperature 
experiments one has to prevent defects from phase ordering by, e.g., stabilizing 
them in a rotating cryostat [26], or trapping them in a superconducting ring [27]. 

KZM for classical phase transitions was confirmed by numerical simulations of 
the time-dependent Ginzburg-Landau model [28] and successfully tested by ex- 
periments in a wide range of condensed matter systems including liquid crystals 
[29], superfluid ^He [26], both high- Tc [30] and low- Tc [27] superconductors, non- 
equilibrium convection systems [31], and Bose-Einstein condensation driven by 
evaporative cooling [32]. With the exception of superfluid ^He - where the early 
detection of defect formation [22] was subsequently attributed to vorticity intro- 
duced by stirring [23], and the situation remains unclear - experimental results 
are consistent with KZM, but more quantitative experimental tests are needed to 
verify e.g. the KZ scaling of defect density with transition rate. 

KZM is a universal theory of the dynamics of continuous symmetry-breaking 
phase transitions whose applications range from the low temperature Bose-Einstein 
condensation to the ultra high temperature transitions in the grand unified theories 
of high energy physics. However, the zero temperature quantum limit remained un- 
explored until recently and quantum phase transitions are in many respects qual- 
itatively different from transitions at finite temperature. Most importantly time 
evolution is unitary, so there is no damping, and there are no thermal fluctuations 
to initialise the symmetry breaking. 

In the next Section, following Ref. [36], we rederive KZM at zero temperature 
where the scalings turn out to be formally the same but the underlying physics is 
different. 
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2.3. KZM in a quantum phase transition 

Near an isolated quantum critical point between two gapped phases both the re- 
action time r and the correlation (or healing) length ^ diverge as 

r~|e|— , (7) 
^~H-^ (8) 

where e is a dimensionless parameter which measures distance from the critical 
point. For instance, when the transition is driven by varying a parameter g in the 
Hamiltonian across a critical point at a finite Qc, then 

6 = ^ . (9) 

9c 

The reaction time determines how fast the system can react to external pertur- 
bations and the healing length sets the scale on which the order parameter heals 
in space, i.e., returns to its ground state value. In a quantum phase transition the 
reaction time is set by an inverse of a gap A between the ground state and the first 
relevant excited state 

r ~ A-i , (10) 

because this is the shortest time scale on which the ground state can adjust adia- 
batically to a varying e. Near the critical point the gap vanishes as 

A - \e\''' (11) 

implying that the evolution across the critical point cannot be adiabatic. 

We consider a second order quantum phase transition that is crossed at a finite 
rate set by the quench timescale tq: 

e(i) = (12) 

In general, e{t) does not need to be a linear function of t, but here we assume the 
generic case when e(i) can be linearised near the critical point as 

e{t) = |(0) t + 0{t') . (13) 

The quench time in Eq. (12) can be identified as tq = |^(0)|~^. A more general, 
but not quite generic, non-linear quench which cannot be linearised near e = is 
considered in Section 2.5. 

Initially, at t — >■ — oo, the system is prepared in the ground state. As long as the 
reaction time in Eq. (7) is fast enough or, equivalently, the gap in Eq. (11) is large 
enough, the state of the system follows its adiabatic ground state. The adiabaticity 
fails near an instant t = —t when the transition rate |e/e| = l/|t| equals the gap 
A ~ lel*^^ = \t/TQ\''^ or 

(14) 
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From this time on long wavelength modes cease to be adiabatic. In a first approx- 
imation, after —t the evolution becomes impulse, i.e., the state effectively freezes 
out between —t and +t when the reaction time r is too slow for the system to 
follow the evolving parameter e(t). The adiabatic evolution of the state restarts 
again at +t when the gap becomes less than the transition rate, see Fig. 1. 
Near the freeze-out time —t, corresponding to 

(15) 

the state is still in the instantaneous ground state with correlation length 

I ~ e-^ ~ r^' . (16) 

In the adiabatic-impulse-adiabatic approximation, this state does not change be- 
tween —t and t. When the adiabatic evolution restarts near t, corresponding to — e, 
the ground state frozen at e becomes an initial excited state for the last adiabatic 
stage of the evolution. 

In this approximation, the impulse stage of the linear quench in Eq. (12) is 
equivalent to a sudden quench from e to — e, where the initial ground state at e 
becomes an initial excited state for the following adiabatic evolution after — e. 

Note that when the quench time tq is large, then e in Eq. (15) is small and the 
linearisation in Eq. (13) is self-consistent because all the non-trivial KZM physics 
happens in the narrow interval between e and — e which is very close to the critical 
point. 

In the adiabatic-impulse-adiabatic approximation the state after the transition 
at — e is equal to the ground state of the system at -|-e. In the adiabatic limit of 
large tq this ground state is very close to the critical point so, by the usual scaling 
hypothesis of the rcnormalisation group theory, expectation value of an operator 
O in this state is proportional to a power of the diverging correlation length ^ in 
the ground state at e. Moreover, since ^ itself scales with the transition time, see 
Eq. (16), the expectation value (O) also scales with a power of the transition time 
TQ. For example, if the phase after the transition admits quasiparticle excitations, 
then their density scales like 

nex ~ , (17) 

where d is the number of space dimensions. The same scale ^ determines a range 
of correlations and other physical quantities. 

In the same way as ^ is the universal scale of length, the freeze-out time t in Eq. 
(14) is the universal scale of time. For instance, the density of excitations during 
the transition depends on time as 

nex(t) ^ F{t/i) , (18) 

where F is a non-universal scaling function. Examples can be found in Refs. [33-35] 
and in Section 2.19. 

KZM predicts a correlation length ^ and a characteristic timescale t. These are 
the basic quantities from which one can derive many other physical observables sim- 
ply by dimensional analysis. This conjecture is confirmed by a series of examples, 
notably the Ising model, where quantities such as density of excitations, excitation 
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energy, correlation functions (often equivalent by a Fourier transform to the mo- 
mentum distribution that is usually measured in the time of flight experiments), 
residual magnetization, entropy of entanglement, generalized entanglement. Berry 
phase, penetration depth in a transition in space, or threshold velocity in an in- 
homogeneous transition are all constructed out of ^ and/or t. These are the basic 
building blocks from which one can construct other physical quantities. 

In this Section wc presented KZM in its most "canonical" form. However, the 
following Sections provide examples how the adiabatic-impulse-adiabatic approxi- 
mation, which is the essence of KZM, can be bend and twisted in many different 
ways as required by actual application. 



2.4. KZM: adiabatic transition in a finite system 

The discussion in the previous Section assumes an infinite system where the corre- 
lation length (16) can diverge to infinity as the transition time becomes infinitely 
slow. However, in a large but finite system of linear size L the divergence must be 
terminated when ^ ~ L, or equivalently 

TQ ~ . (19) 

For slower transitions the scaling (16) breaks down and the transition becomes 
adiabatic. This adiabatic regime is a consequence of a non-zero gap Ac at the 

critical point of a finite system. When the transition is slow enough, then the finite 
gap suppresses any excitation exponentially. This effect was discussed in Ref. [36] 
and it is illustrated by an exact solution in the quantum Ising chain, see Section 
2.13.2 and Refs. [37, 38]. A similar crossover to an adiabatic regime, which can be 
characterized by a steeper scaling of Uex with tq instead of an exponential decay, 
was discussed in Ref. [39] in case of a "half-quench" considered in Section 2.6 below. 



2.5. KZM in a non-linear quench 

The adiabatic-impulse-adiabatic approximation of Section 2.3 can be also applied 
to the non-linear quench considered in Refs. [40, 41], where 



e{t) 



sign(t) 



(20) 



near the critical point e = 0. This function cannot be linearised as in Eq. (13) but, 
nevertheless, essentially the same argument can be applied here as in the linear 
case. 

Indeed, the transition rate |e/e| = r/|t| equals the gap A ~ \e\^^ at e ~ 
{f /'TqY^^^^^^^'^ corresponding to the KZ correlation length and density of exci- 
tations 

I ~ , nex ~ (21) 

respectively. These equations reduce to the corresponding Eqs. (16,17) for the linear 
quench (12) where r = 1. 

Equation (21) shows that the linear quench is not the best choice if we want to 
minimize riex for a given transition time rg, but it is better to take a non-linear 

exponent r S> 1/vz such that nex ~ '^q^^^ ■ This density is less than the density 
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after the linear quench by a factor 

nex(r = 1) ~ 



(22) 



which tends to zero for large tq. 

The non-linear quench with a sufficiently large r S> l/iyz is a good choice if 
we want to improve adiabaticity of the transition. However, the non-linear quench 
requires better experimental control over a system than the linear quench. Not only 
e(t) has to be made non-linear, but it has to be non-linear precisely at the critical 
point. If not, i.e., when the critical point is misplaced by 6e, 



e{t) = 5e — sign(i) 



t 



(23) 



then the quench (23) can be linearised as in Eq. (13) and the linearisation is self- 
consistent for large enough tq when e in Eq. (15) is much less than Se. However, as 
we will see in the next Section and especially in Section 2.20 this technical problem 
does not arise in some "half-quenches" that are robust enough to begin precisely 
at a critical point and go into a gapped phase. 



2.6. KZM: from a critical point into a gapped phase 

In the previous Sections we considered adiabatic passage across an isolated quan- 
tum critical point separating two gapped phases, where KZM was essentially the 
adiabatic-impulsc-adiabatic approximation. Here we consider a transition from a 
critical point into a gapped phase and argue that it can be described by an impulse- 
adiabatic approximation. 

As above, the distance from the critical point is measured by the dimensionless 
parameter 

- (4)' >- ' 

which is in general a non-linear function of time running from to oo. 
The initial ground state at the critical e = has infinite correlation length while 

the ground state at a finite e > has a finite correlation length ^ ~ €~'^. On one 
hand, in first approximation the two ground states look different on length scales 
longer than ^, but they appear the same on distances less than the correlation 
length. On the other hand, the critical ground state is a non-stationary excited 
state with a finite density of excitations with respect to the Hamiltonian at the 
finite e > 0. This excited state and the ground state at e > appear different when 
we look on a scale much longer than the average distance between the excitations, 
but they appear the same when we focus on a shorter scale. Consequently, we can 
identify the density of excitations as Hex — 

The evolution with the non-linear ramp (24) is initially non-adiabatic. It becomes 
adiabatic near t when the quench rate |e/e| = r/\t\ equals the gap A ~ e^'^, or 
equivalently at 

e = e{i) ~ . (25) 

In the crudest impulse-adiabatic approximation, the initial critical ground state 



(24) 
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at e = remains the state of the system until e when the evolution becomes 
adiabatic. Since there is a finite correlation length ^ ~ e"" in the ground state of 
the Hamiltonian at e, the density of excitations is 

nex - . (26) 

It does not change in the following adiabatic evolution after e. The density (26) 
was obtained for the first time in Refs. [39, 42] within the adiabatic perturbation 
theory reviewed in Section 2.12. 

The scaling exponent in Eq. (26) is the same as in the corresponding Eq. (17) 
for a full passage across a critical point, but there is in general a difference in 
numerical pre-factors omitted at the excuse of using We can expect the tt-cx 

in Eq. (26) to be less than that in Eq. (17) because in the "half-quench" starting 
from the critical point the impulse stage is shorter than in the full passage across 
this point. 

However, a quench from a critical point into a gapped phase is not always a 
mere half-quench, as illustrated in Section 2.20 by the process of loading a one- 
dimensional Bose gas (critical Luttinger liquid) into an optical lattice potential. An 
infinitesimally weak potential induces an energy gap proportional to the potential 
strength. This is a realistic experiment, see Ref. [8], where a transition into a gapped 
phase begins exactly at a critical point. 

We will make one more comment on the half-quench after Eq. (49) below where 
we discuss "one-half" of the Landau-Zener transition. In the next Section we de- 
scribe a more unusual application of KZM to a phase transition that takes place 
in space rather than happens in time. 



2.7. KZM in space 

References [44-48] considered a symmetry breaking "phase transition in space" 
where, instead of being time-dependent, a local parameter e(r) is a time- 
independent function of spatial coordinates r . This is a generic scenario in ultracold 
atom gases confined in magnetic/optical traps, where the trapping potential results 
in an inhomogeneous density of atoms p(r) and, in general, the critical point gc 
depends on the local atomic density. Thus even a perfectly uniform parameter g 
translates into an inhomogeneous 

_ 9-gc[pif^] . . 

^ ^ " 9c[p{rl] ' ^''^ 

with the critical point on the surface where e(r) = 0. The part of the atomic cloud 
where e(r) > is in a different phase than the part where e(r) < 0. The phase 
transition between the two phases coexisting in a trap takes place near the critical 
surface e{r) = 0. A classic example of such a phase coexistence is shown in Fig. 2. 
Another plausible experimental scenario can be found in Section 2.21 and Fig. 34. 

In order to apply KZM in this situation, we proceed in a similar way as in Eq. 
(13) and linearise 

€{x) ^ a {x — Xc) , (28) 

near the critical point at Xc where e(xc) = 0. Here the gradient a of e(r) is along 
the X-axis. The system is in the broken symmetry phase where x < Xc and in the 
symmetric phase where x > Xc- In the first "local density approximation" (LDA), 
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Figure 2. Density profiles in the ground state of the one-dimensional Bose-Hubbard model (229) at differ- 
ent repulsion strengths U [43]. The model describes bosonic atoms in an optical lattice potential subject 
to additional harmonic trap confinement, compare Fig. 22 for a typical experimental set up. Here p(i) is 
average number of atoms at site i. The dashed lines are to draw attention to density plateaux at the integer 
value of p = 1 manifesting areas occupied by incompressible Mott insulator phases. These Mott phases 
coexist with superfiuid phases at non-integer p. There are no sharp boundaries between the "integer" Mott 
and "non-integer" superfiuid phases. The phase boundaries are rounded off on a finite length scale ^. This 
scale is long enough to round off plateaux expected at p = 2. (Figure from Ref. [43]b) 



we would expect that the order parameter behaves as if the system were locahy 
uniform, i.e., it is non-zero only for x < Xc and tends to zero as {x^ — x)^ when 
X — >• x~ with the critical exponent /3, see Fig. 3 . 

However, the LDA is not consistent with the divergence of the healing length ^ ~ 
|e|~'^ near the critical point. The diverging ^ is the shortest scale of length on which 
the order parameter can adjust to (or heal with) the varying €{x). Consequently, 
when approaching Xc from the broken symmetry side, the local approximation 
{xc — x)^ must break down when the local correlation length ^ ~ [ci{xc — 
equals the distance (xc — x) remaining to the critical point. Solving this equality 
with respect to ^ we obtain 



a 1+" 



(29) 



Prom the point {x — Xc) — — ^ the "evolution" of the order parameter in x becomes 
"impulse", i.e, the order parameter does not change until (x — Xc) — +£, in the sym- 
metric phase, where it begins to catch up with the local e(x) > and decays to zero 
on the length scale of Thus the "adiabatic-impulse-adiabatic" approximation in 
space, or the "KZM in space", predicts that the non-zero order parameter pene- 
trates into the symmetric phase to the depth of ^ in Eq. (29). As compared to the 
non-analytic LDA prediction {xc — x)^ , in KZM the order parameter is effectively 
"rounded off" on the length scale of ^, see Fig. 3. 

We expect that other physical quantities, which are normally singular or discon- 
tinuous at the critical point, are also "rounded off" on the scale of ^. In particular, 
the energy gap that in a homogeneous system vanishes like A ~ when ^ di- 
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Figure 3. Phase transition in space from a broken symmetry phase to a symmetric phase. Order parameter 
is shown as a function of position x — Xc with respect to a critical point Xc where t{xc) = 0. In the local 
density approximation (LDA) the order parameter is zero in the symmetric phase and follows (xc — x)^ 
with the critical exponent f} in the symmetry broken phase. It is non-analytic at the critical point x^- In 
contrast, the exact order parameter in the ground state of the system is rounded off on the length scale ^ 
in Eq. (29). A non-zero order parameter penetrates into the symmetric phase to a depth ~ ^. 

verges near the critical point, here should remain finite and scale as ^ 

A ~ ~ . (30) 

This finite gap sharply contrasts with the LDA, where one might expect gapless 
quasiparticle excitations localized near the critical point at Xc- 

Two examples supporting KZM in space are described in Sections 2.13.11 and 
2.21 below. Encouraged by these examples, in the next Section we consider an 
inhomogeneous phase transition that, in a sense, takes place in both space and 
time. 



2.8. KZM and dynamics of an inhomogeneous phase transition 

As pointed out already in the finite temperature context [49, 50], in a realistic 
experiment it is difficult to make e exactly homogeneous throughout a system. For 
instance, in the classic superfluid ^He experiments [26] a phase transition was 
forced by neutron irradiation of helium 3. Heat released in each fusion event, 
n + ^He — ?■ ^He, created a bubble of normal fiuid above the superfiuid crit- 
ical temperature Tc- Thanks to quasiparticle diffusion, the bubble was expanding 
and cooling with a local temperature T{t,r) = exp(— r^/2L't)/(27rDt)^/^, where r 
is a distance from the centre of the bubble and D is a diffusion constant. Since this 
T{t, r) is hottest in the centre, the transition back to the superfluid phase, driven 



^A non-linear generalization of Eq. (28), e{x) = sign(x — Xc) \a{x — Xc)\^ , results in a critical point 
rounded off on a length scale ^ ~ q,— i^r/{l-t-i^r) ^ finite energy gap scaling as A ~ a'^'^r/d+vr) ^ 
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by an inhomogeneous parameter 

*r) = 3^4z^, (31) 

proceeded from the outer to the central part of the bubble with a critical front 
rdt), where e(i,rc) = 0, shrinking with a finite velocity v = dr^/dt < 0. 

A similar scenario is likely in the ultracold atom gases in magnetic/optical traps. 
The trapping potential results in an inhomogeneous density of atoms p{f) and, in 
general, a critical point gc depends on atomic density p. Thus even a transition 
driven by a perfectly uniform g(t) is effectively inhomogeneous, 

(32) 

9Mr)] 

with the surface of critical front, where e{t, r) = 0, moving with a finite velocity. 

According to KZM, in a homogeneous symmetry breaking transition, a state after 
the transition is a mosaic of finite ordered domains of average size Within each 
finite domain the orientation of the order parameter is constant, but uncorrelated to 
orientations in other domains. In contrast, in an inhomogeneous symmetry breaking 
transition [49] the parts of the system that cross the critical point earlier may be 
able to communicate their choice of orientation of the order parameter to the parts 
that cross the transition later and bias them to make the same choice. Consequently, 
the final state may be correlated at a range longer than ^, or even end up being 
a ground state with long range order. In other words, the final density of excited 
quasiparticles may be lower than the KZ estimate in Eq. (17) or even zero. 

Prom the point of view of testing KZM, this inhomogeneous scenario, when 
relevant, may sound like a negative result because an imperfect inhomogeneous 
transition suppresses system excitation by KZM. However, from the point of view 
of adiabatic quantum computation or adiabatic quantum state preparation, it is the 
KZM itself that is a negative result: no matter how slow the homogeneous transition 
is, there is a finite density of excitations (17) in the final state which decays only 
as a power of the transition time tq. Prom this perspective, the inhomogeneous 
transition may be the way to suppress KZ excitations and prepare the desired final 
ground state adiabatically. 

In order to estimate if and when the inhomogcneity is actually relevant, we 
linearise the parameter e{t, x) in both t and x near the critical front where e(t, x) = 
0, as 

e{t, x) a (x — vt) , (33) 

in a similar way as in Eq. (12). Here a is an inhomogeneity of the quench and v is 
velocity of the critical front. When observed locally at a fixed x, the inhomogeneous 
transition in Eq. (33) looks like the homogeneous quench in Eq. (12) with 

rg = —■ (34) 

av 



The part of the system where x < vt, or equivalently e{t, x) < 0, is already in the 
broken symmetry phase. The orientation of the order parameter chosen in this part 
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can be communicated across the critical point not faster than a threshold velocity 

■0 ~ ^ . (35) 

When V ^ V the communication is too slow for the inhomogeneity to be relevant, 
but when v <^v we can expect the final state to be less excited than predicted by 
KZM. 

Given the relation (34), the condition (35) can be rewritten either as 

■0 a i+- , (37) 
or as a relation between the slope and the threshold transition time, 

TQ a 1+" . (38) 



The last relation means that, for a given inhomogeneity a, the transition is effec- 
tively homogeneous when tq <^ tq, but the inhomogeneity becomes relevant when 
the transition is slow enough and tq^ tq. In the homogeneous limit of a — > the 
threshold transition time tq — > oo. 

An example of inhomogcncous transition is solved in detail in Section 2.13.12 
and in Rcfs. [44, 45, 48, 51, 52]. In the next Section we return to the basics and 
review the Landau-Zener model. 



2.9. The Landau-Zener (LZ) model 

In this Section we briefly review the Landau-Zener (LZ) model [17] because in a 
number of integrable models KZM can be derived exactly by mapping a model 
to a set of independent LZ anti-crossings, see Section 2.11 and the Sections that 
follow. After this brief review in Section 2.10 wc will attempt to invert the relation 
between the LZ model and KZM and estimate LZ excitation probabilities in the 
adiabatic-impulse-adiabatic approximation essential for KZM [53, 54]. 

In dimensionless variables, the LZ model is defined by a two-level time-dependent 
Hamiltonian 

H = ') (39) 



2 V 1 -<t) 
where 

e{t) = ±- (40) 

and TQ is a transition time. The Hamiltonian passes through an anti-crossing at 
e = 0. At any fixed e, it has two instantaneous eigenstates: the ground state | | 
(e)) and the excited state j t (e))- In the time-independent basis |1),|2) of the 
Hamiltonian (39), the instantaneous eigenstates are 

|t(e)) = |l)cos^ + |2)sin^ , \ i (e)) = - |1) sin ^ + |2) cos ^ , (41) 
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Figure 4. In panel A, the instantaneous energy spectrum of the two-level system in Eq. (39) parametrised 
by the time t. In panel B, the adiabatic {t < —i or t > i) and impulse regimes (—i <t<i)in the two-level 
system dynamics. Note the similarity of panel B to Figure 1. (Figure from Ref. [54]) 



where cos e = e/^TT? and sm9 = Xj^YT^. Their eigenenergies axe A/2 and 
—IS. 1 2 respectively, where 



A 



(42) 



is an instantaneous energy gap. It is minimal at the anti-crossing centre at e = 0. 
A good example of the Landau- Zener two level system is the Feshbach resonance 
in Figure 5, see also Section 2.23. 

In the LZ problem the system is prepared in the ground state at an initial time 
ti and we want to know the probability P that it is excited at a final time t/. A 
general solution to this problem was obtained in Ref. [17] and then analysed in 
Refs. [54, 56]: 



m)) = 



\l)[2idt + ^^+\2) 



[a L>_i_j^g/4 (iz) + b Z>_i_j^g/4 {-iz)] (43) 



where z = ;^/=e Dm{s) is a Weber function [57], and a, b are constants to be 
determined by initial conditions. 

The excitation probability P can be obtained in a closed form in a number of 
special cases: 

(i) In the textbook case of ti — oo and tf ^ oo we obtain the exponential 
LZ formula 



(44) 



(ii) When the evolution begins at = in the ground state at the anti-crossing 
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— closed channel 

— open channel 




Figure 5. The plot shows a scattering potential for two atoms as a function of their separation r. The 
scattering potential depends on internal states of the colliding atoms. For a small scattering energy E the 
atoms can be either in an open (thick solid line) or closed (thin solid line) channel. When one of the bound 
states in the closed channel with energy Em is close to the scattering energy E Ri 0, then we have the 
Feshbach resonance [55] . An external magnetic field coupled to magnetic moments of the colliding atoms 
can shift the energy in the closed channel moving the energy Em either up or down with respect to the 
open channel. The field can be used to tune the scattering atoms close to the Feshbach resonance [55] . In 
Section 2.23 we consider scattering of two fermions, corresponding to fermionic annihilation operators c-f 
and cj^, near resonance with a bound molecular state Em, correspoding to a bosonic annihilation operator 
b. A coupling between the two fermions and the molecular state is decribed by a Landau- Zener Hamiltonian 

H = — ^eb^b+ jt (^c|c^ + c|c^^ + § (^^^c^c^ + "^l"^!^) ■ Here e is a distance from the resonance proportional 
to a difference between the actual and resonant values of the magnetic field. The state with two fermions 
corresponds to the state |1) and the state with a molecule to the state |2) in Eq. (39). In a scattering of two 
fermions far from the resonance a weakly occupied molecular state can be elliminated and in a second order 

perturbative expansion one obtains an effective Hamiltonian for fermions H^g = i ^c|c^c|c4_^ . Depending 

on the sign of e this is additional repulsion or attraction, so the magnetic field can be used to change the 
strength and even sign of the interaction between the fermions. This effective interaction strength diverges 
at € = 0, but this is where the occupation of the molecular state is large so it cannot be elliminated and 
has to be included explicitly in the calculations. 



centre and runs to tf ^ oo, then according to Refs. [54, 56] 



2sinh(^) _ 



TTTq/ 



TTTQ 




1 

4 



(45) 



where T{x) is the gamma function [57] and the last approximate form is 
the asymptote when tq S> 1. Note that, unlike in the full LZ transition (i), 
here the probability does not decay exponentially with tq, but only as Tq^. 

(iii) When the evolution runs from tj — )■ — oo to = 0, then the excitation 
probability P is the same as in case (ii), see Ref. [54, 56]. 

(iv) When the LZ problem is not symmetric with respect to the anti-crossing 
and 



eit) = 



— , when t <0 

TQ 

, when t>0 



(46) 



with S ^ 1, then the excitation probability after a passage from ti — >■ — oo 
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to tf ^ +00 is [54, 56] 



P=l - - sinh — ^ e i — 



2 V 4 

1/1-5 



2 



' Tq' , (47) 

where the last asymptote is accurate when rg S> 1. Note that when S ^ oo 
the second half of the transition becomes adiabatic, the excitation proba- 
bility does not change in this adiabatic part of the evolution, and the final 
P becomes the same as in case (iii) where the transition ends at t/ = 0. 

In the adiabatic limit of rg ^ 1, the excitation probability P is exponentially 
small in the standard LZ model (i), but in the non-symmetric cases (ii,iii,iv) there 
is much slower power-law decay P ~ ^q'^- these cases an instantaneous rate of 
the transition de/dt is a discontinuous function of time. For instance, in case (ii) 

de _ JO,whent<0, 

Tt = \:^,whent>0. ^^^^ 

The discontinuity at t = translates into a fat high frequency tail in the Fourier 
transform of e(i) which helps to excite the system despite its finite energy gap. 
More generally, as discussed in Refs. [17, 39], any non-analytic behaviour of e(t) 
results in a power law decay of the excitation probability. A discontinuity in the 
r-th derivative of e{t) leads to ^ 

P ^ Tq'^ (49) 

for large enough tq. 

Case (ii) above is similar to the half-qucnch in Section 2.6. This connection will 
become more transparent in Section 2.11, where KZM will be represented as a 
series of independent LZ transitions. The non-linear parameter e{t) in Eq. (20) has 
a discontinuous r-th time derivative at t = 0. Since on a microscopic level KZM 
excitations originate from LZ transitions, then one might expect their density to 
decay as ricx ~ Tq^^ for large enough tq. However, in the thermodynamic limit 
the density of excitations riex is dominated by contributions from low frequency 
modes, with wavelength longer than ^, whose excitation probabilities are close to 
1 and for whom the asymptotic tail in Eq. (49) does not apply. 

These comments complete our brief review of the LZ model. In Section 2.11 below 
the LZ model is applied to derive KZM in a class of exactly solvable models, but in 



-'-Indeed, in the basis of instantaneous eigenstatcs (41) the state is \tp{t)) = o| t (f))e 2*/-oo )1 _|_ 

j3\ ]. {e))e~^^'' -^-"^ '^['!{* ^jjg adiabatic limit we have |a| ^ l,/3 ~ 1 and, to leading order in the 

small a, we obtain a perturbative excitation probability 



/ dt'^{U^\i)e 
J_oo dt de 



2 



The discontinuity in Eq. (48) yields P ~ Tq for rg 2> 1 • A more general non-linear 

, , _ f , when t <0 
'^^ ' ~ \ W'^qY .when t > 

has a discontinuous r-th derivative at t = and we obtain P ^ Tq'^^ for large tq. 
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the next Section we will attempt to rederive and reinterpret cases (i-iv) using the 
adiabatic-impulse-adiabatic approximation essential for KZM, see Refs. [53, 54]. 



2.10. The Landau- Zener model in the adiabatic-impulse approximation 

A unitary evolution with the time-dcpcndcnt Hamiltonian (39) is adiabatic when 
t —> ^oo and the gap in Eq. (42) is large enough, but it may be not adiabatic 
near the anti-crossing centre at t = where the gap is minimal. In the adiabatic- 
impulse-adiabatic approximation, the evolution is adiabatic before —t and after 
+t, but it is impulse between —t and +t, see panel (b) in Figure 4. The crossover 
time t between the adiabatic and impulse stages is the time when the transition 
rate measured by |e/e| = l/\t\ equals the instantaneous energy gap in Eq. (42): 

where a 2± 1 is an adjustable parameter, compare Fig. 4. The solution is 



" - A - J_ 
^ ~ TQ ~ V2 

With this estimate at hand, we can reiterate the cases (i-iv) listed in Section 2.9: 

(i) The evolution (39) starts from the instantaneous ground state at ti <C — t 
and ends at t/ S> t. In the adiabatic stage before —t the state follows the 
instantaneous ground state. Then in the impulse stage between —t and +t 
the state does not change and remains equal to the instantaneous ground 
state I \. (— e)) at — t. Thus in the adiabatic-impulse-adiabatic approxima- 
tion, the excitation probability at t is 



'1 + 



- 1 



(51) 



Pai = Kt(6)U(-e))P = (52) 

and it does not change in the following adiabatic evolution after t. Its 
expansion up to second power of tq, 

Pai = 1 - aTQ + + 0{4) , (53) 

matches the corresponding expansion of the exact exponent in Eq. (44) 
when we set a = 7r/2 ~ 1. 
(ii) The evolution (39) begins at = in the instantaneous ground state 
I 4 (0)) and ends at S> £ Since this state does not change in the initial 
impulse stage of the evolution until t, the excitation probability at t is 

Pai = |(t(e)U(0))P = ^ (l - -^=L=) , (54) 

and it does not change in the following adiabatic evolution after t. Its 
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expansion in powers of tq is 

^Ai = ^ - Iv^+liargf/' + Oir'^/') , (55) 
while the corresponding expansion of the exact solution (45) is 

Given that 2 — 41n2/7r = 1.1 1, the two expansions are in good agreement 
when we choose a = tt/A ~ 1. 

(iii) The evolution (39) begins at ti <^ —t and ends at t/ = 0. It is adiabatic 
before —t and impulse from —t to 0. In the impulse stage the state remains 
the instantaneous ground state | J, (— e)) at —t. Similarly as for the exact 
solution, the final excitation probability Pai = Kt (0)1 4- (— e))P is the same 
as in case (ii). 

(iv) The evolution (39) begins in the ground state at ti <^ —t and ends at 
tf S> ts- The transition is not symmetric: when t > the transition time is 
tqS instead of tq for t < 0. Consequently, the crossover from the adiabatic 
to impulse stage takes place near —t as before, but the second crossover 
from the impulse to adiabatic stage is at ts given by Eq. (51), but with tq 
replaced by tqS. Similarly as in case (i), the excitation probability is 

PAi = \{Ue5)\i{-m' 

= 1 - i (l + v^)' arQ + 1(1 + 6){1 + V~S)\aTQf + 0(r|) ,(57) 

while the corresponding expansion of the exact solution (47) is 

P = l-l(l + V5yTQ + ^{l + 6){l + V6fT'Q + Oi4). (58) 

The two expansions match when we set a = 7r/2 2± 1. 

In all cases, after fitting only one adjustable parameter a ~ 1, we obtain accurate 
leading and next-to-leading order terms of the expansion in powers of tq, so the 
adiabatic-impulse-adiabatic approximation is accurate when tq ^ 1 and the ex- 
citation probability P is substantial. Since the KZM excitation density nex in Eq. 
(17) is dominated by contributions from long wavelength modes, whose excitation 
probability is close to 1, the adiabatic-impulse-adiabatic approximation can accu- 
rately predict nex- The approximation is not accurate for long wavelength modes, 
but their contribution to the total nex is negligible. 

In view of the above conclusion, it is not quite surprising that in the next Section 
the LZ model provides exact solutions which are consistent with KZM. 

2.11. KZM as a set of independent Landau-Zener transitions 

The KZM argument in Section 2.3 was confirmed by exact solutions in a number 
of intcgrablc models, sec e.g. Rcf. [35, 37, 38, 58, 59] and the following Sections. 
There is, however, an assumption in Section 2.3 that the transition passes through 
an isolated quantum critical point between two gapped phases. As demonstrated 
by exact solutions in a number of integrable models, KZM does require careful 
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generalisation for a transition across a multicritical point [60, 61], along a gap- 
less line [62], or across a gapless phase [63]. The generalisation was obtained in 
integrable spin models, such as the XY spin chain or the two-dimensional Kitaev 
model [64], all of which can be mapped to non-interacting fermions by a Jordan- 
Wigner transformation [18]. A transition in a translationally invariant system of 
non-interacting fermions can be mapped to a set of independent LZ anti-crossings. 
The general LZ argument summarized in this Section was gradually developed in 
Refs. [35, 37, 53, 58-63, 65-68]. 

We assume a general local Hamiltonian quadratic in fermionic cre- 
ation/annihilation operators and translationally invariant on a d-dimensional lat- 
tice. Its quasimomentum representation is 

where A; is a quasimomentum running over the first Brillouin zone, is fermionic 
annihilation operator, and is a 2 x 2 Hermitian matrix. The Hamiltonian can 
be diagonalized by a Bogoliubov transformation 

where 7^ annihilates a fermionic Bogoliubov quasiparticle, and (ug, v^) is an eigen- 
mode of stationary Bogoliubov-de Gennes equations 



(^) 













(61) 



with positive eigenfrequency w^. The diagonalized Hamiltonian is H = 
^^cj^ (^l^^lj: — 5) • Its ground state is a Bogoliubov vacuum |0) annihilated by all 
7^. This ground state is an initial state for a dynamical quantum phase transition. 

A time-dependent problem can be solved in the Heisenberg picture, where the 
state remains the initial Bogoliubov vacuum |0), but the operators eg evolve in time. 
Their evolution can be conveniently described by time-dependent modes («^, Vj^) 
in Eq. (60). Indeed, given that 7^ remains time- independent in the Heisenberg 
picture, a Heisenberg equation i^cj^ = [cj^,H(t)] is equivalent to time-dependent 
Bogoliubov-de Gennes equations 

Thus the problem was separated into a set of independent two-level systems enu- 
merated by k. 

For each k, an initial state [u^(— 00), 00)] is the positivc-cj^ cigcnmode of Eq. 
(61) with the initial H^(—oo). This positive eigenmode represents the initial vac- 
uum ground state for the quasiparticle 7^. What we want to know is a probability 
Pj^ that after the dynamical transition a quasiparticle 7^ of the final Hamiltonian 
H{oo) is excited. Since this excited state is represented by the negative eigenmode 
of Eq. (61) with eigenfrequency —ujj:, then is equal to the probability that the 
solution of Eq. (62) ends in the negative eigenmode of the final iyg(oo) or, sim- 
ply, the initially excited two-level system becomes deexcited. This is essentially the 
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same question as in the Landau-Zener (LZ) model in Section 2.9 so in order to 
answer this question we map Eq. (62) to the LZ model in Eq. (39). 

Near a quantum phase transition driven by a parameter e the two-level Hamil- 
tonian can be linearised, or is linear, in the small e: 

H^{t) = e{t) a{k) + a'{k) . (63) 

Here iy{k) = ^^aj(/c)cr*, ct* with i = x,y,z are Pauli matrices, cr'{k) = ^io!,i{k)(T^ , 
Ui and a'j are model-dependent functions, and e{t) = t/rq is the dimensionless 
parameter in Eq. (12) quenched from t — oo to t ^ +oo. Notice that in this 
general Hamiltonian we ignored a c-number term proportional to the 2x2 identity 
matrix bccaTisc it would contribute to the global phase only, but see Ref. [69]. The 
instantaneous eigenvalues of the Hamiltonian if^ are iw^, where 



oj^ = ^J{ed + a'f . (64) 

The question is what is the probability for the initially excited state to become 
deexcited to the ground state — wg? 

To answer this question, we map the general Eqs. (62,63) to the LZ model (39). 
The main obstacle is that, unlike in the LZ model, the "spin matrices" a and a' 
in the Hamiltonian (63) are not necessarily pointing in orthogonal directions. This 
minor problem can be easily fixed by orthonormalisation ^. The orthonormalised 
Hamiltonian is 

H^{t) = \a{k) e{t) + h{k)\ a{k) + A{k) a±{k) (65) 
with instantaneous quasiparticle spectrum 

ooj; = ^[e a{k) + b{k)^^ + A^k) . (66) 

We can eliminate b(k) by shifting the time variable in e{t) = t/rq to t' = t + 
TQ h{k)/a{k): 

t' - - - - 

HAf) = — a{k) a{k) + A{k) a±{k) . (67) 

TQ 

Up to an unimportant rotation of "spin" quantization axes, this is the Landau- 
Zener Hamiltonian (39) and we can use the LZ formula for the (de)excitation 
probability 



-"-Indeed, -wc rewrite the matrix o-(k) = a(k)a-(k), -where a = \/W is the length of vector a, and ct(A;) = 
a(fe)o"' is a normalized Pauli spin matrix pointing along the unit vector a(k) = d{k)/a{k). Next, we 
decompose cr' {k) = b{k)a{k) + A(k)a±(k), where b{k) = d(k) d'(k) is the component of vector a'(fc) along 
the unit vector a(fe), a±{k) = a'{k) — b{k)a{k) is the component of a'(fe) orthogonal to a{k), A(fe) = 

y 3±{k)'^ is the length of the orthogonal component, and &± = "^ff (t' is a normalized Pauli spin 

matrix pointing perpendicular to &{k). 
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k 

Figure 6. The minimal gap function A(A;) = |2sin(fe)| and excitation probabiUty pj. = exp 

for different quench times tq in the one-dimensional transverse field quantum Ising chain in Section 2.13. 
The gap function is zero at the Fermi point kp = 0. With increasing tq the excitation probability becomes 

a Gaussian p^j exp(— 27rTQfc^) localised around the Fermi point. Its integral ricx = / ^Pk = Tq 

is the density of excited quasiparticles. 



2 



Consequently, the density of quasiparticles is an average of over the first Brillouin 
zone 

1 V- f d'^k 

= NX.Pk ^ J J^Pk^ (69) 

where is a number of lattice sites and the integral becomes accurate in the 
thermodynamic limit N ^ oo. 

In order to evaluate the integral in Eq. (69), we note that in the LZ Hamiltonian 
(67) the minimal instantaneous gap at t' = is twice A{k). It vanishes on a "Fermi 
manifold" (point /line/surface) in quasimomentum space where 

A(fc) = . (70) 

Consequently, Eq. (68) implies that = 1 only on the Fermi manifold and, in the 
adiabatic limit rg — )■ oo, the integrand in Eq. (69) remains non-negligible only 
very close to the Fermi manifold, see the examples in Figures 6 and 7. In this limit 
the integral (69) can be easily done in a few special cases: 

(i) When the quench runs through an isolated critical point between two 
gapped phases, as in Section 2.3, then there is an isolated Fermi point 
kp where A{kF) = 0. When the Fermi point is isotropic in /c-space, then 
the gap vanishes like 

A2(fc) ~ \k-kF\'^ (71) 
near kp, and the integral in Eq. (69) yields 



(72) 
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(a) (b) 

Figure 7. In the two-dimensional Kitaev model on a honeycomb lattice in 
Section 2.15 the minimal gap function is A(fc) = 2[Ji sin(A;Mi) — J2 sin(A;M2)], 

where Mi = j + and M2 = — f i are spanning vectors of triangular 
reciprocal lattice. Here we assume J^; = Jy = 1 so that the Fermi line is 
ky = 0. The 3D plots show excitation probabilities pg for tq = 1 (left panel) 
and TQ = 16 (right panel). With increasing tq the excitation probability 
localizes on the Fermi line. 



when Tq is large enough for to be locahsed close enough to kp. Examples 
of case (i) are the transverse field quantum Ising chain in Section 2.13 and 
the quench across a multicritical point of the transverse field XY chain in 
Section 2.14. 

(ii) More generally, when an isolated Fermi point in d dimensions is not isotropic 
and has different critical exponents in different directions, then the density 
of excitations scales as 



nex ~ Tq^'+- (73) 

for large enough tq. Here the exponents z, v and apply in respectively 

m and (d—m) directions. This is the case of the semi-Dirac points in optical 
lattices considered in Ref. [67]. 

(iii) When the Fermi manifold is an m-dimensional surface and the minimal 
gap depends on the distance from the Fermi surface as A'^{k) ~ \k — kp]^^ , 
where kp is the point on the Fermi surface closest to k, then the integration 
in Eq. (69) yields 

^ex ~ Tq^ (74) 

when Tq is large enough. An interesting example is the 2D Kitaev model 
in Section 2.15 and Ref. [63], see Figure 7. 

(iv) When, as in case (iii), the Fermi manifold is an m-dimensional surface with 
A^{k) ~ \k — kp\^^ near a generic point kp on the surface, but there is 
an isolated "dominant" Fermi point kp where /^{k) ~ \k — kp\^^ with 

< ZA; then the integral in Eq. (69) is dominated by a neighbourhood of 
the dominant Fermi point and 

nex ~ Tq'^''^ (75) 

when Tq is large enough. An example can be found in Ref. [35]. 

Case (i) is the problem treated in Section 2.3 by the standard KZM, while cases 
ii-iv) are generalisations that go beyond that simple argument. However, even in 
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case (i) we obtain the density of excitations (72) with an exponent d/zA which 
in principle can be different than the exponent diy/{l + uz) in the KZM equation 
(2.3), unless 

ZA = — • (76) 

Here za on the left hand side characterizes the minimal gap function A(/c) while 
the critical exponents z, v on the right hand side characterize the critical point. 

To see if the relation (76) holds, we return to the Hamiltonian (65) with the 
quasiparticle spectrum (66). Without loss of generality, we can assume that h{kp) = 
at the Fermi point ^ and its asymptote near kp is 

b^ik) - \k-kF\^' . (77) 

On one hand, the critical exponents z and u are defined by the asymptotes near the 
Fermi point: ~ {k — kpl^ at e = and lo^ ~ \e\'^^ for small e, see Ref. [15]. On the 

other hand, the actual asymptotes of the spectrum (66) are ~ — 2™"^(^'"^^) 
at e = and ojj^^ ~ |e|^ for small e. Comparing the definitions with the actual 
asymptotes we obtain z = ^mm{zb, za) and uz = 1. The equality (76) holds if 

ZA < Zb (78) 

or, equivalently, near the Fermi point the function b^{k) is negligible as compared 
to A^(fc). The inequality is the condition for the KZM argument in Section 2.3 to 
give the same scaling as the exact solution in case (i). 

This condition may be not quite surprising because when we neglect 6^(fc), then 
for each k the instantaneous gap in Eq. (66) is minimised at the critical point e = 0. 
Since the critical point is the anti-crossing center for each k, the dynamical excita- 
tion of each two-level system takes place near the critical point, so the excitation 
probabilities must be determined by the critical exponents that characterize the 
quasiparticle spectrum near this critical point. Given the accumulated evidence, 
the condition (78) seems to be satisfied by most quenches across an isolated quan- 
tum critical point, but it is not fulfilled by quenches across a multicritical point 
considered in Refs. [60, 61] and Section 2.14. 

For instance, in the quench considered in Section 2.14 one finds za = 6, = 4 
and z = l/i' = 2, the inequality (78) is not satisfied, and the exact scaling ricx ~ 
Tq^^^ in Eq. (72) is different than the Uex ~ '^Q^^'^ predicted by Eq. (17). In this 

example 6^(fe) is not negligible, so it shifts the anti-crossings away from the critical 
point, and the excitation probabilities are dominated by non-critical energy 
modes away from e = 0. The shifts, and the anomalous exponent 1/6, originate 
in the fact that the dynamical excitation process takes place asymmetrically with 
respect to the multicritical point. As a consequence, dynamical scaling may require 
introducing new non-static exponents ^ 



-'This can be always achieved by a fc-independent shift e' = e + b{kp)/a{kp). 

^Indeed, the effective Hamiltonian (67) has instantaneous spectrum ix>L = ^ (e')^ + A^, where e' = t' /tq. 
Prom this spectrum we formally obtain effective exponents z' = za/2 and v' z' = 1. With these effective 
exponents z' and v' in plax;e of the "canonical" z and v, the KZM Eq. (17) gives the exa<;t scaling in Eq. 
(72). 
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In this Section we obtained exact scalings by mapping a class of integrable 
fermionic systems to the LZ model. In the next Section we rederive the KZM scal- 
ing in Eq.(17) from the adiabatic perturbation theory. This approximate derivation 
does not seem limited to integrable models. 



2.12. KZM from adiabatic perturbation theory 

An alternative derivation of the KZ scaling in Eq. (17) was presented in Refs. 
[39, 70]. In the integrable models equivalent to non-interacting fermions considered 
in the previous Section, this derivation is a perturbative approximation to the 
exact LZ argument, but its validity does not seem limited to integrable models. 
The argument is as follows. 

Let the set of functions 4>p{e) represent an eigenbasis of a Hamiltonian -ff(e) 
depending on a parameter e. A wave function of the system can be expanded in 
the eigenbasis 

^ = T.cip{r)Me) ■ (79) 



Wc assume the linear e(t) = t/rq in Eq. (12) with a large tq — >■ oo. The Schrodinger 
equation is equivalent to 

i^ + ^^(^q{p\-^\'l)=^P<^P ^ (80) 

*3 q 

where u}p{e) is the eigenenergy of the Hamiltonian H{e) corresponding to the eigen- 
state 0p(€) (or \p)). After a unitary transformation eliminating the dynamical phase 

ap{t) = ap{t) e-'^*'^^' '^"l^^*')! = ap(e) e''^'^ '^"^^'^ , (81) 
the Schrodinger equation becomes 

^ = -Y^aq{e){p\^\q) e^^^'*' . (82) 

etc (xc 

? 

Since the system were initially prepared in the ground state |0), then in the adia- 
batic limit the single term g = dominates the sum and the excitation probability 
is given by 



/oo J 
de — 10) fJ-TQ r de'[>^„{e')-^^{e')] 



(83) 



The derivative of the ground state ^|0) is finite at a continuous phase transition. 
We assume a uniform d-dimensional system with a single (relevant) branch of 

excitations which can be labelled by (quasi-)momentum k. The branch has a dis- 
persion relation with a gap A(e). The gap is non-zero everywhere except the critical 
point at e = 0, where it vanishes like A ~ \€\'^^ and the excitation A; = is gapless. 
Conservation of momentum implies that only pairs of excitations with opposite 
momenta {k, —k) can be excited by the time-dependent e. In this framework, Eq. 
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(83) gives the following density of excitations 



29 



lip. 



(27r)' 



Pk 



d 



(84) 
(85) 



Here \k,—k) is an excited state with one pair of excitations of opposite momenta 
±k. 

A general scaling argument implies that 



u^{e)-uo{e) = AF{A/k') = F{\er /k') , 



(86) 



where k = \k\ and F (or F) is a non- universal function with a universal tail 
F{x) 1/x for large x. Substituting Eq. (86) to Eq. (85) we obtain 



Pk 



dx {k,-k\^\0) e' 
dx 



(87) 



where x = k "^l^t and x' = k ^^'^e'. Another scaling argument 



\zu—l 



-fcl^lO) = sign(6)^^G(H-A^) 



transforms Eq. (87) into 

f°° J • / \ I \zu—l i\ \zv\ ''■('''Qk^^^ f^dx Ix'l"" Fdx'l"") 

w / ax sign (x) I x I G{\x\ ' ' ' i ^ 



, (89) 



which depends on k only through the dimensionless combination rqk « . This 
means that there is a characteristic wave vector 



k 2± Tr,^^ 



(90) 



corresponding to the KZ wavelength ^ ~ fc ^ in Eq. (16). The k marks a crossover 
between a regime of large k ^ k^ where the adiabatic approximation is self- 
consistent and Eq. (89) accurately predicts small p^, and a regime of small k <^ k 
where we do not expect Eq. (89) to be accurate. It is expected to be an overestimate 
for fermionic excitations and an underestimate for bosons. 

Since k is the only momentum scale in Eq. (84), then for dimensional reasons 
the integration yields 



no 



t—d ^ 1+21^ 



(91) 



up to possible logarithmic corrections [70]. This is again the KZ equation (17). 

This Section completes our general discussion of KZM. The following Sections 
provide specific examples in a number of model systems. 
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2.13. Quantum Ising chain: transition across an isolated critical point 

Much of our understanding of quantum phase transitions is based on the proto- 
typical integrable quantum Ising chain [15] 

N 

H = - E (5 < + (92) 

n=l 

with periodic boundary conditions un+i = o"i . In the thermodynamic hmit 
N ^ oo the model has two critical points at = =tl between a ferromagnetic 
phase when 1^1 < 1 and two paramagnetic phases when \g\ > 1. Their critical 
exponents are z = u = 1. Today the quantum Ising chain is not only an exactly 
solvable toy model, but it is also becoming a subject of quantum simulation in 
experiments: the NMR simulation in Ref. [11] and the ion trap experiment in Ref. 
[9], see Figs. 8 and 9 respectively. 
Here, like in Refs. [36-38, 58], we consider a linear ramp 

git) = -- (93) 

with t running from — oo to across the critical point at = 1- The quench 
begins at 5 — > oo in the ground state | ttt • • • t) with all spins polarized along 
the z-axis. At the final g = there arc two degenerate ferromagnetic ground states 
with all spins pointing either left or right along the x-axis: | — )•— )•— )• • • • — )■) or 
I ••• ■h-). In an adiabatic classical transition from the paramagnetic to 

ferromagnetic phase, the system would choose one of the two ferromagnetic states. 
In the analogous quantum case, any superposition of these two states is also a 'legal' 
ground state providing it is consistent with other quantum numbers conserved by 
the transition from the initial paramagnetic state. 

However, when — t- 00, then energy gap at g = 1 tends to zero (quantum 
version of the critical slowing down) and it is impossible to pass the critical point 
at a finite speed without exciting the system. As a result, the system ends in a 
quantum superposition of states like 

with finite domains of spins pointing left or right and separated by kinks where 
the polarisation of spins changes its orientation. Average size of the domains or, 
equivalently, average density of kinks depends on the transition rate. When the 
transition is slow, then the domain size is large, but when it is very fast, then 
orientation of individual spins can become random, uncorrelated with their nearest 
neighbours. 
When we define a dimensionless parameter 

e = l^ = g-l, (95) 

9c 

then according to KZM in Section 2.3 the evolution is adiabatic before e ~ ^/^^ 
and after — e, and impulse near the critical point. The correlation length in the 
adiabatic ground state at e is proportional to 



(96) 
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Figure 8. NMR gate sequence to prepare an effective pure state | in the quantum Ising chain (92) 

by spatial averaging from thermal equilibrium state. (Figure from Ref. [11]) 




Figure 9. Schematic picture of the ion trap experiment in Ref. [9]. By adiabatically increasing effective 
spin-spin interaction an initial paramagnetic state is driven into a superposition of two ferromagnetic states 
with opposite magnetisation. (Figure from Ref. [9]) 



This ground state is (approximately) the initial state for the last adiabatic stage of 
the evolution after — e. This argument shows that when passing across the critical 
point, the state of the system gets imprinted with a finite KZ correlation length 
proportional to ^. In particular, this coherence length determines average density 
of kinks after the transition as 

nex ^ = Tq'/' . (97) 

This is an order of magnitude estimate with a pre-factor ~ 1. 

In the following we derive an exact solution [37, 38, 58] not only to confirm 
the scaling exponent in Eq. (97) and provide an unknown numerical pre-factor, 
but also to see how much the state after the transition resembles the adiabatic 
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ground state at e. In the adiabatic limit tq — )■ oo this near-critical ground state 
has a divergent correlation length ^ which, by the standard scaling hypothesis of 
the renormalisation group, determines all physical observables according to their 
scaling dimensions. 

2.13.1. Landau-Zener argument 

Here we assume that N is even for convenience. After the non-local Jordan- 
Wigner transformation [18], 

< = (4 + Cn) n (1 - 24cm) , (98) 

Tn<n 

al = i (4 - Cn) n (1 - 2cLc„) , (99) 

m<n 

< = 1 - 2cicn , (100) 

introducing fermionic operators which satisfy |c^,cJi| = Smn and {c^,Cn} = 
|c]ruc|i| = the Hamiltonian (92) becomes [71] 

= P+ P+ + P- H- P~ . (101) 

Above 

P± = i[l±P] (102) 
are projectors on subspaces with even (+) and odd (— ) parity 

N N 

p = n < = n (i - 24cn) (103) 



n=l n=l 



and 

N 

= 5Z (254cn - 4cn+i - - C„+lCn - 4c^l - (104) 

n=l 

are corresponding reduced Hamiltonians. The c^'s in H~ satisfy periodic bound- 
ary condition c^v+i = ci, but the c„'s in obey cn+i = — ci - what we call 
"antiper iodic" boundary conditions. 

The parity P of the number of c-quasiparticles is a good quantum number and 
the ground state has even parity for any g ^ 0- Assuming that a quench begins in 
the ground state, we can confine to the subspace of even parity. is diagonalised 
by a Fourier transform followed by a Bogoliubov transformation [71]. The Fourier 
transform consistent with the antiperiodic boundary condition is 



-in /A 
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where the pseudomomentum k takes "half-integer" values 

127r iV - 1 27r 

The Hamiltonian (104) becomes 

H'^ = 5Z l^ic/ - cos k]clck + sin k cj^cLj^. + c-kCk 



33 



(106) 



-g] . (107) 



Diagonalisation of iJ+ is completed by a Bogoliubov transformation 



Ck = Uklk + V-kl-k ' 



(108) 



provided that Bogoliubov modes {uk,Vk) are eigenstates of the stationary 
Bogoliubov-de Gennes equations 



CO 



2{g — cosk) 2 sin A; 
2sinfc —2{g — cosk) 



Uk 
Vk 



(109) 



There are two eigenstates for each k with energies lj = iwfc, where 



Uk = 2-v/ {g — cos ky + sin k . 



(110) 



The positive energy eigenstate 



iuk,Vk) 



Ok . Ok 

COS — , sm — 
2 ' 2 



(111) 



where tan^^ = smk/{g — cos A;), defines the quasiparticle operator 



,t 



Ik = U^Ck + V_kC'_^. , 



(112) 



and the negative energy eigenstate (uj^ ,v^) = {—Vk,Uk) defines 7^ = (u^^ )*Ck + 
{vZ^T^Lk — ~l-k- After the Bogoliubov transformation, the Hamiltonian 



2 ^k ^« Ok^k - 7fc ^Ik ) is equivalent to 



(113) 



This is a simple-looking sum of quasiparticles with half-integer pseudomomcnta. 
However, thanks to the projection P"*" in Eq. (101) only states with even 

numbers of quasiparticles belong to the spectrum of H. 

In the linear quench (93), the system is initially in its ground state at large 
initial value of _g ^ 1. hwt as g is ramped down to zero, the system gets excited 
from its instantaneous ground state and, in general, its final state at t = has 
finite number of kinks. Comparing the Ising Hamiltonian (92) at 5 = with the 
Bogoliubov Hamiltonian (113) at 5 = we obtain a simple expression for the 
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operator of the number of kinks 



1 ^ 



n=l k 

The number of kinks is equal to the number of quasiparticles excited at = 0. The 
excitation probability 

Pk = iihk) (115) 

in the final state at ^ = can be found with the time-dependent Bogohubov 
method. 

The initial ground state at 5 ^ oo is a Bogoliubov vacuum |0) annihilated by all 
quasiparticlc operators 7^ determined in Eq. (112) by the positive-oj eigenmodes 
[ukiVk) ~ (1)0) of Eq. (109). As g{t) is ramped down, the quantum state gets 
excited from the instantaneous ground state. We work in the Heisenberg picture 
where the state remains the initial Bogoliubov vacuum and the operators do 
not depend on time, but the operators are time-dependent 



Uk{t)lk + v-k{t)llk ' (116) 



with the initial condition [uk{— 00), Vk{— 00)] = (1,0). They satisfy Heisenberg 
equations i-^Ck = [ck,H^] equivalent to a time-dependent version of the 
Bogoliubov-de Gennes equations (109): 



d 

dt V Vk 



Uk 



2[g{t) — cos A;] 2 sin A; 
2 sin k —2[g{t) — cos k\ 



Uk 
Vk 



(117) 



These equations are an example of the general equations (62). 

Comparing equations (62,65) and (117) with g{t) = 1 -|- e{t), we can identify 
a{k) = 2, b{k) = 2(cosA;- 1), 

A(A;) = 2 sin A; , (118) 

(j(A;) = — cr^, and d'±{k) = a^. There are two Fermi points kp = 0,it where A^kp) = 
which correspond to the critical points g'c = 1) — 1 respectively. Modes near kp = 
do not get excited in our quench from y — )• 00 through i^c = 1 to = 0, so we can 
focus on Af = only. Expansion in k — kp like in Eqs. (71,77) yields za = 2 and 
Zfe = 4 and the LZ argument in Section 2.11 case (i) yields 

nex ^ Tq^/- = Tq"^ (119) 

for TQ » 1. This is the same as the KZ prediction (97) because the condition 
< Z}y in Eq. (78) is satisfied. 
What is more, with the minimal gap function A « 2|A;| near kp = Q we have 

« e-^-^Q'^^ , (120) 



localised near A = when tq S> 1, see Figure 6, and the integral in Eq. (69) yields 
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the density of excitations 



nex = -A=- 'Tn''^ . (121) 



The numerical pre- factor is = 0.113. For comparison, in Ref. [70] the 

perturbative Eqs. (84,85) were used to estimate the same density of defects as 

— 1/2 

nex ~ 0.124 Tq . The scaling is the same as in the exact solution (121) but, as 
a result of the perturbative approximation in Eq. (85), the numerical pre- factor is 
overestimated as expected for fermionic quasiparticles, see Section 2.12. 

2.13.2. Adiabatic transition in a finite chain 

No matter how long tq is, a quench across a critical point in an infinite system 
is not adiabatic because the energy gap at the critical point is zero. However, in a 
finite system the gap is finite and we can expect the quench to become adiabatic for 
a sufficiently long tq, see Section 2.4. For instance, in the periodic Ising chain the 
(relevant) energy gap at gc = 1 for excitation of the lowest two energy quasiparticles 
with k = ±Tr/N is w+^/^v + '^-■k/n = 47r/A/" and decays like 1/N. Since pk is a 
probability to excite a pair of quasiparticles {+k, — k) and different pairs evolve 
independently, a probability for the state to follow the adiabatic ground state 
without any quasiparticles is a product 

Pgs = Hil-Pk) (122) 

fe>0 

over the positive "half-integer" k = {m + 1/2)2'k/N in Eq. (106). Well on 
the adiabatic side p-^/n ~ exp (— 27r^-^) ^ 1 is exponentially small and 

P(i+2m)7r/w ~ {P'k/n)'^^^'^^^^ ^ Ptz/n foi" > 0, SO we can safely approximate 

^'gs = 1 - P^/N + O (pIi^) « 1 - exp (-2vr'^^) . (123) 
This Pgs ~ 1 and the quench in a finite chain is adiabatic when 

TQ » ^ , (124) 

i.e., when even the lowest energy pair of quasiparticles with k = ±7r/N is not likely 
to get excited. Reading this inequality from right to left, the size iV of a defect-free 
chain grows like Tq^. This is consistent with Eqs. (96) and (121). 

2.13.3. Exact solution of the time- dependent Bogoliubov-de Gennes equations 

If we want to go beyond the simple result for Uex in Eq. (121), then we need an 
exact solution of Eq. (117). When A; > a new time variable 

T = 4rQsinA; ( — -l-cosA; ) (125) 



brings Eqs. (117) to the canonical LZ form (39) 



dT\Vk) 2 



-RkT 1 

1 RkT 



Vk 



(126) 
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with a transition rate = l/Argsin^ k. A general solution of the Landau-Zener 
equations (126), see Section 2.9, is 

Vkir) = aD^s-ii-iz) + bD^s-i{iz) , 

Ukir) = (^TRk - 2i^^ Vkir) , (127) 

with arbitrary complex parameters a,b. Here Dm{x) is a Weber function, s = 
and z = ■y/RkTc''^^^. The parameters a,b are fixed by the initial conditions 
Uk{—cc) = 1 and Vk{—oo) = 0. Using the asymptotes of the Weber functions when 
r ^ — oo, we obtain 

g-7r/8Rfc 

« = 0, \b\' = ^^- (128) 

At the end of the quench at t = when r = 2rQsin(2A;), the argument of the 
Weber function is iz = ^/Rtre^'^^^ = 2y^e*'^/^ cos(fc)sign(fc). When tq 3> 1 the 
modulus of this argument is large for most k, except near k = and we can 
again use asymptotes of Weber functions to obtain products 

\Vkf = 1 - \Ukf , 



- sin A; + sign(/c) e^''^^'^'' Vl - e-2'^^«'=' e'"^" 



UkVk 2 

<^fe = ^ + 2rQ - (2 - ln4)TQA;2 + A;Vq Iutq - arg [r (l + iTqk'^)] . (129) 

These products depend on k and tq through two dimensionless combinations: Tqk'^, 
which implies the usual KZ correlation length ^ = ^Jtq, and fc^rg In tq which 
implies a second scale of length \/tq Inrg. The final quantum state sX g = cannot 
be characterised by a unique scale of length. Physically, this reflects a combination 
of two processes: KZM that sets up the initial post-transition state of the system 
at — e, and the subsequent adiabatic evolution with quantum dephasing. 

2.13.4- Entropy of a block of spins after a quench 
The Von Neumann entropy of a block of L spins, 

S{L) = - Tt PL log2PL , (130) 

measures entanglement between the block and the rest of the chain. Above pL is 
reduced density matrix of the subsystem of L spins. In recent years this entropy 
was studied extensively in ground states of quantum critical systems [72-79] . At a 
quantum critical point, the entropy diverges like log L for large L with a pre-factor 
determined by the central charge of an effective conformal field theory [72, 76, 80]. 
In the quantum Ising model at the critical = 1 

SGS(L) ~ ilogsL (131) 
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for large L. Slightly away from the critical point, the entropy saturates at a finite 
asymptotic value [76] 



(132) 



when the block size L exceeds the finite correlation length ^ in the ground state of 
the system. 

In a dynamical quantum phase transition the quantum state of the system devel- 
ops a finite correlation length ^ = ^/tq. If this dynamical correlation length were 
the only relevant scale of length, then one could expect the entropy of entangle- 
ment after a dynamical transition to be given by Eq. (132) with ^ simply replaced 
by ^. However, as we could see in Eq. (129), there are two scales of length, and - 
strictly speaking - there is no reason to expect that either of them alone is relevant 
in general. This is why we cannot rely on scaling arguments and the entropy has 
to be calculated "from scratch" . 

We proceed in a similar way as in Ref. [72, 73, 75, 77, 79] and define a correlator 
matrix for the block of L spins 



n 



a, /3t 
/3,1-a 



(133) 



where a and /? are L x L matrices of quadratic correlators 



1 1 e"^^ 



(134) 



and 



sign(m — n) 



^\,\m-n\ 



2iTQ- 



(135) 
(136) 



where we used Eqs. (129). Here ^ = ^/tq is the KZ correlation length and 

I = In TQ (137) 

is a dephasing length ^. When the linear quench is extended to (7 — > — oo, as in 
Section 2.13.10, the dephasing length continues to grow like l{t) = 4t. 



^Indeed, for large tq, the non-ground-state part of I3m,,n, i-c the second term in the square bracket in 
Eq. (136), is negligible on distances \rn — n| <C ( as a result of dephasing the integral over k in Eq. (135). 
Its integrand has the phase ipj. in Eq. (129) whose fast rotation with k is driven mainly by the term 
k^TQ liiTQ. This term introduces a length scale Zq = y'TQ In tq and, if it were the only length scale, the 
integral would be negligible or "dephased" for |m — n| ^ lo- However, Eq. (129) depends on k also through 
the combination rgfc^ introducing ^ as the second scale. Due to conspira<;y between the two scales, the 
integral is actually dephased on distances up to the net dephasing length I = Iq/^ in Eq. (137) 
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We note that a and P are Toeplitz matrices with constant diagonals. The ex- 
pectation values (...) are taken in the Bogoliubov vacuum state. Since this state 
is Gaussian, all higher order correlators can be expressed by the matrices a and 
/3 - they provide complete characterization of the quantum state after the dynam- 
ical transition. Since the matrices depend on both scales ^ and both scales are 
necessary to characterize the Gaussian state. 

As noted in Ref. [72, 73, 75, 77, 79], the block entropy can be calculated as 

S{L,tq) = - Trp log2P = - IVn log^U. (138) 

In this calculation we use Eq. (128) and Eqs. (134,135). The calculation involves nu- 
merical evaluation of the integrals in Eqs. (134,135) and numerical diagonalisation 
of the matrix 11. Results are shown in Panel A of Fig. 10. The entropy grows with 
the block size L and saturates at a finite value Soo{tq) for large enough L. In Panel 
B, the asymptotic entropy is fitted with the (linear) function Soo{tq) = A+BIutq. 
On one hand, KZM suggests a simple replacement of the ground state correlation 
length ^ in Eq. (132) by the dynamical KZ length ^ = implying the asymptotic 
value 

Sooirq) ^logal ^ Inrg = 0.120 Ihtq . (139) 

On the other hand, the best fit gives B = 0.128 ± 0.004 and A = 1.80 ± 0.05. The 
best B is in reasonably good agreement with the expected value of 0.120. 

In Panels C and D of the same figure, the entropy S(L,tq) is rescaled by its 
asymptotic value -S'oo(tq) Ri A + Blnrg. After this transformation we can better 
focus on how the entropy depends on the block size L. A simple hypothesis would 
be that the entropy depends on ^ and saturates when L ^. To check if this 
is true, in Panel C the block size L is rescaled by ^ = ^^tq and this rescaling 
brings the plots close to overlap, although they do not overlap as well as one might 
have hoped. By contrast, as shown in Panel D, rescaling the block size L by the 
dephasing length I makes the multiple plots collapse. Thus we can conclude that 
the entropy saturates at 

5oo (tq) ~ ^ logsl when L » I (140) 

i.e. the entropy of a large block of spins is determined by the KZ dynamical cor- 
relation length ^ = ^/TQ, but the entropy saturates when the block size is greater 
than the dephasing length I = ^/tq In tq . A similar conclusion was also obtained in 
the numerical simulations in Ref. [81], where the non-integrable case with non-zero 
longitudinal magnetic field was considered as well. 

According to KZM the correlation length ^ is determined in the impulse stage 
when the system is crossing the critical point, while the second scale I builds up 
by dephasing the excited state in the following adiabatic stage. This scenario is 
confirmed by the block entropy at the moment when the system is crossing the 
critical point at = 1, see Figure 11. The entropy saturates at 

Socirq) ^ ^ logsl when L > |. (141) 

Near the critical point, when the scale I set up by dephasing just begins to build 
up, ^ is still the only relevant scale of length. The result (141) confirms the KZM re- 
placement rule ^ — )• ^ in Eq. (132). As predicted by the adiabatic-impulse-adiabatic 
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Figure 10. Panel A shows entropy of a block of L spins after the dynamical phase transition as a function 
of the block size L. The multiple plots correspond to different values of the quench time tq. For all rg, 
the entropy grows with the block size L and saturates at a finite value Sco{tq) for large enough L. In 
Panel B, this asymptotic value of entropy is fitted with the function Soo{tq) = A-\- Blnrq. The best fit 
has B = 0.128 it 0.004 and A = 1.80 it 0.05. This B is in reasonably good agreement with the expected 
value of B = = 0.120. In Panels C and D, entropy S{L, tq) is rescaled by the best fit to its asymptotic 
value 5oo (tq ) = A + BItitq. With this rescaling one can focus on how the entropy depends on the block 
size L. In Panel C, the block size is rescaled by ^ = y'rQ but the rescaled plots do not overlap. However, 
as shown in Panel D, rescaling the block size L by the second scale I = yPrq In tq makes the six plots 
collapse. (Figure from Ref. [38]) 



approximation in Section 2.3, the entropy is the same as if the state were the ground 
state at e. 

2.13.5. Impurity of a block of spins after a quench 

Since a fully analytic calculation of entropy does not seem possible, it is worth- 
while to calculate another more easily tractable entanglement-related quantity. For 
example, an "impurity" of the correlator matrix H in Eq. (133) 

/(n) = Tr n (1 - n) (142) 

is zero only when the L spins are in a pure state i.e. when all eigenvalues of 11 are 
either or 1. It is maximal when all the eigenvalues are ^, or when the state is 
most entangled. Thanks to its simple quadratic form, the impurity can be calculated 
relatively easily. 

Simple calculation using the block structure of 11 in Eq. (133) and the Toeplitz 
property of the block matrices a and /3 leads to 

I{L) =2|Lao- (L-|j|)(aj2 + |/3,f) I , (143) 
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Figure 11. Panel A shows entropy of a block of L spins during the dynamical phase transition at the 
critical point g = 1 as a function of the block size L. The multiple plots correspond to different values of 
the quench time tq. For all the quench times, the entropy grows with the block size L and saturates at 
a finite value 5oo (tq) for large enough L. In Panel B this asymptotic value of entropy is fitted with the 
function Soo{tq) = A + Blnrg. The best fit has B = 0.126 ± 0.005 and A = 0.80 ± 0.03. This B is in 
reasonably good agreement with the expected value of iJ = = 0.120. In Panels C and D, the entropy 
S{L,tq) is rescaled by the best fit to its asymptotic value Soo{tq) ~ A + Blnrq. With this rescaling 
one can focus on how the entropy depends on the block size L. In Panel D, rescaling the block size L by 
^ = y/TQ makes the six plots collapse. (Figure from Ref. [38]) 

where aj = aj^o and 13 j = Further calculation as in Ref. [38] shows that at the 



final 5 = the impurity saturates at 



— J- when L S> / , 



(144) 



in analogy to the entropy in Eq. (140). 

It is interesting to compare the dynamical impurity (144) with the impurity in 
the ground state near the critical point. A simple calculation gives the asymptote 
of impurity at the critical point 



InL 



(145) 



when InL ^ 1. Near the critical point, the asymptote is valid when the block is 
much shorter than the correlation length, L ^ ^, and at larger L the impurity 
saturates at 



rGS 



In^ 



(146) 



Again, in the same way as for the block entropy, the simple KZM replacement rule 
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— 7- ^ applied to Eq. (146) gives the correct asymptotic value of the dynamical 
impurity in Eq. (144). As predicted by the adiabatic-impulse-adiabatic approxi- 
mation in Section 2.3, the dynamical impurity is the same as the impurity in the 
ground state at e. 

2.13.6. Correlation functions after a quench 

Correlation functions are of fundamental interest in the theory of phase transi- 
tions because they provide direct manifestation of their universal properties and 
are in general directly accessible experimentally. Here we present results for spin- 
spin correlation functions after a dynamical quantum phase transition, see Ref. 
[38, 58]. 

To begin with, we observe that for symmetry reasons the magnetization (a^) = 0, 
but the transverse magnetization 

«) = (1 - 2cicn) = 2ao - 1 ~ 2'K^/2¥^ ^^^^^ 

when TQ 3> 1. This is what remains of the initial magnetization (a^) = 1 in the 
initial paramagnetic ground state at 3 — > oo. As expected, when the linear quench 
is slow, then the final magnetization decays towards (cr^) = characteristic of the 
ferromagnetic ground state at y = 0. 
Final transverse spin-spin correlation function at 5 = is 

Cf = - «)«+H> = ■i(\flR?-WR?) « 



■kR-' 
' 2 1^ 



(148) 



when R> 1 and Inrg ^ 1. This correlation function depends on both ^ and the 
dephasing length / in Eq. (137), but its long range tail 

Cf ~ e"^ (149) 

decays on the scale I. 
In contrast, the ferromagnetic spin-spin correlation function 

Cl^ = {«+r) - K)K+r) = K<+r) (150) 

cannot be evaluated so easily. As is well known, in the ground state can be writ- 
ten as a determinant of an i? x Tocplitz matrix whose asymptote for large R can 
be obtained with the Szego limit theorem [82]. Unfortunately, in time-dependent 
problems the correlation function is not a determinant in general. However, below 
we circumvent this problem in an interesting range of parameters. 
Using the Jordan- Wigner transformation, can be expressed as 

C^^ = {boaibia2 . . . 6^-10^) . (151) 

Here a„ and 6„ are Major ana fermions defined as Un = (cl + c„) and bn = — Cn. 
Here a„ is hermitian and 6„ is anti-hermitian. Using Eq. (134) and Eq. (135) we 
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Figure 12. Ferromagnetic spin correlation in Eq. (153). The oscillatory term means that consecutive kinks 
are anti-correlated - they keep more or less the same distance ~ ^ from each other forming something similar 
to a ...-kink-antikink-kink-antikink-... regular lattice with a lattice constant ~ ^. However, fluctuations in 
the length of bonds in this lattice are comparable to the average distance itself giving the exponential 
decay of the correlator on the same scale of ~ ^. 



obtain 



(152) 



Using Wick theorem, the average in Eq. (151) is a determinant of a matrix when 
{amCLn) = and (bmbn) = for m 7^ n, or equivalently when Im/3m-n = for 
m ^ n. Inspection of Eq. (136) shows that Im/3m-n ~ when \m — n\ <^ I. 
Consequently, when R <C / then we can neglect all Im/3n2_j^ assuming that {cLm(^n) — 
and {bmbn) = for m ^ n. The correlation function is a determinant of the 
Toeplitz matrix [{bman+i)]m n=i R within the quasiparticle horizon where R <^l. 
Asymptotic behaviour of this Toeplitz determinant can be obtained [38] using 
standard methods in Refs. [82] and [58, 83]: 



oc 



exp 



-0.174 5 



cos 




(153) 



when 1 <C i? <C / i.e. up to the dephasing length /, see Fig. 12. 

The final ferromagnetic correlation function (153) exhibits decaying oscillatory 
behaviour on length scales much less than the scale I, but both the wavelength of 
these oscillations and their exponentially decaying envelope are determined by ^. 
As discussed in a similar situation by Cherng and Levitov [58], this oscillatory be- 
haviour means that consecutive kinks keep more or less the same distance ~ ^ from 
each other forming something similar to a ...-kink-antikink-kink-antikink-... lattice 
with a lattice constant ~ ^. However, fluctuations in the length of bonds are com- 
parable to the lattice constant itself and they give the (marginally underdamped) 
exponential decay of the correlator on the same scale of ~ ^. 

The gaps in analytic knowledge of the correlation functions were filled by nu- 
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Figure 13. Panel A shows the dynamical transverse correlation Cj^ as a function of magnetic field g in 
the linear quench. For each tq, we show both numerical (dashed) and analytical (solid) result. The plots 
overlap near the critical point at g = 1 but diverge in the ferromagnetic phase when g < 1 indicating a 
breakdown of our numerical simulations in this regime. Panel B shows analytic and numerical results for 
the dynamical transverse correlation function at the moment when the quench crosses the critical point at 
g = 1. The transverse correlators overlap well confirming that our numerical simulations are still accurate 
at the critical point. Finally, in panel C, we show the dynamical ferromagnetic correlation function 

at g = 1 and in panel D, we show the same correlation function after rescaling R/^. The rescaled plots 
overlap quite well supporting the idea that near the critical point the KZ correlation length ^ = y^rg is 
the only relevant scale of length. (Figure from Ref. [38] ) 



merical simulations with the real time TEDB (time evolving decimation block) 
algorithm [84] in Ref. [38]. As illustrated in panel A of Figure 13, the simulations 
were stable enough to cross the critical point and enter the ferromagnetic phase, 
but once in the ferromagnetic phase, the algorithm was breaking down due to the 
quasiparticle horizon effect, see Section 3.2. The numerics can be trusted at = Ij 
but it is not reliable when g < 1. KZM can be tested at the critical point, but one 
cannot reliably follow the dephasing in the ferromagnetic phase. 

Panel B of Figure 13 shows the transverse correlation at gc = 1 for several 
values of tq. For each tq, both the numerical correlator and its analytic counter- 
part from Eq. (148) are plotted and they approximately coincide. Equation (148) 
can be also used to obtain analytically, but with some numerical integration, the 
exponential tail of the transverse correlator when tq 3> 1: 



Czz 
R 



0.44 



exp 



R 

-2.03^ 



(154) 



accurate when -R ^ ^. This tail decays on the KZ correlation length ^ which proves 
to be the relevant scale of length. 
Panel C shows the ferromagnetic spin-spin correlation functions at = 1 for 
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the same values of tq. They arc roughly exponential and their correlation length 
seems to be set by ^ = ^/tq. To verify this scaling hypothesis in panel D the same 

plots as in panel C are shown but with R rescaled by ^. The rescaled plots overlap 
reasonably well confirming the expected ^ = ^/fq scaling. 

The correlation functions at the critical point = 1 are consistent with the 
prediction of the adiabatic-impulse-adiabatic approximation in Section 2.3 that the 
state in the impulse stage is approximately the ground state at e with correlations 
decaying exponentially on the correlation length ^. 

2.13.7. Fidelity between the state after a quench and the final ground state 

A fidelity F = \ {iIj\GS)'\^ between the final state and the desired final ground 
state |GS) is a very sensitive measure of the quality of adiabatic quantum com- 
putation or adiabatic quantum state preparation. In the quantum Ising chain, the 
fidelity can be obtained analytically from the exact solution, sec Rcf. [85]. The 
fidelity F is a probability that not a single pair of quasiparticles with opposite 
quasimomenta {k, —k) is excited after a quench, 

F = Hil-Pk), (155) 

fe>0 

where the product runs over the positive half-integer quasimomenta in Eq. (106). 
It is more convenient to calculate its logarithm 



InF = Vln(l-pfe) « / dkln{l-pk), 
k>0 



(156) 



accurate for ^. If the transition terminates aX g = 0, then for tq S> 1 we have 
Pk = exp(— 27rrQ/s^) and 

Nhp-x fn° ds In Tl — e~*n 
InF = - ^oo,^ ^ = -1.3ne.Ar. (157) 

2 Jo dse ' 

Here ngx = Jl^ ^Pk = ^/'^t^\/'^'^Q is the density of excitations in Eq. (121). Thus, 
as might have been expected, the fidelity decays exponentially with the chain size 
N on the length scale ^ ~ n~J^, 

F = e'^-^ , (158) 

but with an unexpected coefficient 1.3 greater than 1. 

It is instructive to compare the decay of fidelity in Eq. (158) to a simple Poisso- 
nian model where each of the N bonds is either excited (with probability ricx) or 
not excited (with probability 1 — ricx) independently of other bonds. In this model 
the fidelity is a probability that none of the N independent bonds is excited: 

i^Poisson = ( 1 - nex )'^ ~ e'^^^ (159) 

when Hex ^ 1- Comparing Eqs. (158) and (159) we find that the decay in Eq. 
(158) is faster than in the model of independent bonds. The faster decay means 
anti-bunching correlations between the kinks randomly scattered along the spin 
chain in the final state. The same anti-bunching correlations are also responsible 
for the oscillations in the ferromagnetic correlation function in Eq. (153) and Fig. 
12. 
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The anti-bunching in Eq. (158) contrasts with the bunching apparent in fidehty 
between the ferromagnetic ground state at < \g\ < 1, which has finite density of 
kinks riex, and the fully polarized ferromagnetic ground state at ^ = 0: 

i^ferro = e"'"-"^ , (160) 

sec Ref. [85] . Here the coefficient is | meaning that the relevant density of uncorre- 
lated excitations is one-half of the density of kinks. The kinks in the ferromagnetic 
ground state are bound into kink-antikink pairs randomly scattered along the chain. 
The bound pairs do not destroy the ferromagnetic long range order. This is most 
apparent when g <^ 1 is perturbatively small and the ground state is a fully polar- 
ized ferromagnet but with an admixture of a small density of reversed spins. Each 
reversed spin is a tightly bound pair of kink and antikink which, however, does not 
affect the long range ferromagnetic order. 

2.13.8. Geometric phase after a quench 

Connection between the geometric phase [86] and quantum phase transitions has 
been explored recently in Refs. [87]. The geometric phase can be used as a tool 
to probe quantum phase transitions in many body systems. It is a witness of a 
singular point in the energy spectrum arising in non-trivial geometric evolutions. 
In Ref. [88] the geometric phase of the final state at = was calculated as 

r = 2'kN (riex — 1) — Stt -I- TT ^cos TTHex + sin -kN cot • (161) 

It depends on the density of excited kinks ricx in Eq. (121). 

2.13.9. Generalized entanglement in the final state 

The generalized entanglement of a state ]V') with respect to a Lie algebra h with 
hermitian, orthonormal basis of generators O; enumerated by I = 1, M can be 
characterised by a purity of {ip) relative to h which is given by 

M 

1=1 

where A/" is a normalisation factor chosen so that < < 1. 

A fermionic system, like e.g. the quantum Ising model in the representation 
of Jordan- Wigner fermions in Eq. (104), can be characterised by generalized en- 
tanglement with respect to the algebra h of number preserving bilinear fermionic 
operators spanned by c|cj with 1 < i,j < N. In Ref. [89] it was shown that, near an 
isolated quantum critical point, the generalized entanglement of the ground state 
|GS(€)) scales like 

PU|GS(e))]-P41GS(0))] ~ , (163) 

where ^ is a correlation length in the ground state. The relative entanglement with 
respect to the critical point is inversely proportional to the diverging correlation 
length. 

Moreover, as verified in Ref. [89], in a linear quench across an isolated critical 
point we have 



(164) 
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where \ip{t)) is the state of the system during the transition, |GS[e(i)]) is the in- 
stantaneous ground state, and G is a scaUng function. The relative amount of en- 
tanglement with respect to the instantaneous ground state is inversely proportional 
to the KZ length ^, which diverges in the adiabatic limit, and its timc-dcpcndence 
becomes universal when measured in units of t in Eq. (14), as expected from KZM 
in Section 2.3. 

In particular, the simple replacement rule ^ — t- ^ applied to the generalized 
entanglement in the ground state (163) gives the correct dynamical result (164) 
in accordance with the adiabatic-impulse-adiabatic approximation in Section 2.3 
where the state during the impulse stage is approximately the ground state at e 
with correlation length ^. 

2.13.10. Linear quench to t +oo: quantum dephasing after KZM 

As we could see in Section 2.13.6, near the critical point gc = ^ the state of the 
system is close to the adiabatic ground state at e which can be characterised by a 
single KZ correlation length ^ = ^/fQ, but in the following adiabatic evolution to 
g = quantum dephasing develops the second longer scale I = ^^tq In tq. In this 
Section, as in Ref. [58], we extend the linear quench to t — t- +oo to give the quantum 
dephasing enough time to complete its job, and see if there is any imprint of the 
KZ length ^ left in the final state at t — >■ -|-oo after the dephasing is completed. 

The extended quench crosses both critical points gc = ±1 and the quasiparticle 
excitation probability in Eq. (68) becomes a sum of two Gaussians localised near 
the two Fermi points kp = 0, tt: 

= e-^'rrQsin^fe ^«>>1 ^-2nr^k^ ^ ^-2^Mk-nr _ (jgg) 

Consequently, the average density of excitations in Eq. (69) , 

= ^ ^q'^' ' (166) 

is twice the density at 3 = in Eq. (121). This density of t spins excited 

in the final ground state | X-IX ■ ■ ■ i) at — >■ —00. 

Another consequence of crossing two critical points and Eq. (165) is a simplified 
quadratic correlator in Eq. (134), 

which is zero when m — n is odd. The other correlator Prn,n in Eq. (135) is even 
simpler. 

Indeed, the last adiabatic stage of the evolution, after crossing the second critical 
point at gc = —1, leaves plenty of time for quantum dephasing. When g <^ —1 
the eigenstates of the Bogoliubov-de Gennes Eqs. (109) are {uk,Vk) ~ (1,0) and 
(ufe, Vk) ~ (0, 1) with eigenvalues — and respectively, where u^, ~ 2(5 — cos k). 
Consequently, the adiabatic solution of Eqs. (117) is 

Uk ~ VPke^'i , Vk ^yi - Pk e . (168) 



As t — >■ 00 the solutions {uk,Vk) accumulate phase factors exp^z(2tcos A; — t^/rq) 
which oscillate strongly with k. These oscillations "dephase" the integral in Eq. 
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(135) to zero: 

Pm,n = ^ J[ dk VW^) ^'^-"^-'-''-^ e^M-«) (169) 
when t — > oo. 

However, when t is large but finite, then the dcphasing in Eq. (169) is effective 
only when \ra — n\ <^ 4t, where the 4 is twice the maximal group velocity of 
instantaneous quasiparticles when \g\ 3> 1. This estimate means that the dephasing 
length I in Eq. (137) continues to grow like 

l{t) = At (170) 

in the last adiabatic stage of the evolution after gc = —1. 

Notice that the correlator am,n depends on ^. Since in the absence of any non-zero 
/5m, n the quadratic correlator am,n = {cmcli) contains all the information about the 
final Gaussian state, then we can conclude that ^ remains a permanent imprint of 
KZM which has not been washed out by quantum dephasing. This imprint has 
several interesting manifestations: 

• Using Szego limit theorem, the ferromagnetic correlator can be found as [58] 

CW oc exp (-0.174|) cos - ^^o) (171) 

when TQ » 1 and 1 <C i? ^ /(t). As a consequence of crossing two critical points 
instead of one, this correlator is zero for odd R. When t — >■ oo then l(t) — >■ oo 
and the correlator becomes accurate for all i? ^ 1. 

As we can see, during the linear quench the ferromagnetic correlator evolves 
from the pure exponential decay at g^c = 1, see Figs. 13C and D, to the oscillating 
exponential decay at g ~oo, sec Eq. (171). These oscillations arc gradually 
developed by quantum dephasing during the adiabatic stages of the evolution. 
At 5 = the oscillations are limited to i? <C y/rqlnrQ, but as g ^ — oo they 
gradually extend to all (even) R. The dephasing relaxes the state to a crystal-like 
train of kinks and antikinks. 

• The transverse magnetization is 

«) = 2ao,o - 1 = - 1 + 2nex , (172) 

i.e., the magnetization —1 in the final ground state at y — )• — oo plus twice the 
density of inverted spins in Eq. (166). 

• The transverse correlator is 

= - 4|ao,il| ~ ^ 7r2"|2~ ' ^ ' 

compare Eqs. (148,166,167). This is roughly the Z — >■ oo limit of the same corre- 
lator at g = 0, compare Eq. (148). 

• The block impurity introduced in Section 2.13.5 is 



7(n) = Trn(l-n) = 2TVa(l-a) ^ 2n^^L (174) 
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Figure 14. In A and B, spontaneous ferromagnetic magnetization as a function oin — Uc and the rescaled 
X = y/a{n—nc) respectively. The collapse of the rescaled plots in panel B demonstrates that the penetration 
depth of the ferromagnetic order parameter into the paramagnetic phase is 5n ~ a~^/^ which is the same 
as the ^ predicted in Eq. (29). (Figure from Ref. [45]) 



when L > ^ > 1 but L < l{t). Here we used Eqs. (133,166,169,167). 

• The extensive impurity in Eq. (174) implies that the block entropy is also ex- 
tensive: S{L,tq) = sL for L S> ^. The entropy density s is most conveniently 
evaluated as the entropy of the whole chain S{N, tq) divided by A^, see Ref. [58]: 

-2Tralog2(l-a) 
s = lim 

N^oo N 

2 dk 

^ ~hi2 y 27r [Pfe^°S2PA: + (1 -Pfc)log2(l -Pfc)] ~ 0.11 riex (175) 

accurate when n^x ^ 1- Here we used Eqs. (133,138,167,169) and evaluated the 
trace in the quasimomentum basis taking the continuous k limit for N ^ oo. 

Notice that here the entropy of the whole chain sN appears finite while the 
actual quantum state of the isolated chain is pure. This is an artifact of sending 
/3m, n — )• in Eq. (169) which is accurate only when |m — n| <^ l{t) at a finite t. 
This condition limits validity of the extensive entropy S = O.llricxL to L <^ l(t). 

• Two site entanglement at g ^ — oo was studied in Ref. [90]. Both concurrence 

and negativity are zero when the distance R between the two sites is 
odd. Only entanglement between even-neighbour sites can be non-zero. For large 
Tq, these two measures of entanglement scale as 

^ ^'-1 , ~ , (176) 

when n is fixed, and they decay exponentially with n/^ for large n. What is 
more, the concurrence and negativity are non-zero only when tq > Tq""^ , where 

Tq^^ is a threshold that increases with n. When tq < Tq^ there is no two-site 
entanglement and the entanglement is entirely multipartite. 

2.13.11. Transition in space 

In this Section we support the general discussion in Section 2.7 by a solution in 
the quantum Ising chain. We consider the ground state of an open Ising chain in an 
inhomogeneous transverse field g„ which can be linearised near the critical point 
at Uc where = 1 as 



ffn 1 + a (n - ric) , 



(177) 
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compare Eq. (28). The open chain is in the ferromagnetic phase where n < ric 
and in the paramagnetic phase where n > ric- We want to know if the non-zero 
spontaneous ferromagnetic magnetization = (cr^) in the ferromagnetic phase 
penetrates across the critical point into the paramagnetic phase and what is the 
depth of this penetration? 

The quadratic Hamiltonian iJ+ in Eq. (104) with open boundary condi- 
tions is diagonaUsed to = '^"^'^'^'Tm by a BogoUubov transformation 
c„ = X]m=o(^"™^™ "^mnlm) with m enumerating N eigenmodes of stationary 
Bogoliubov-de Gennes equations 

^mut,m = '^9nU^,m " 2<^i,^ (178) 

with LOm > 0. Here = Unm =t Vnm- We make a long wavelength approximation 
^nq=i m ~ ^n,m T 'Sii''^n,m ^-^d obtain a long wavelength differential equation 

UmU^ = 2a{n - nc)u^ ± "^dnu"^ . (179) 
Its eigenmodes with Um > are 

UJm = VSma ,Um{n) OC V'm-l(a;) +'il)rn{x) ,Vm{n) (X tpm-lix) - ''Pm{x) , (180) 

where 

X = \/a{n — He) (181) 

is a rescaled distance from the critical point, V'm>o(a^) are eigenmodes of a harmonic 
oscillator satisfying ^(—5^ + .T^)'0m(x) = (m + l/2)i/jrn{x) , and ip-i{x) = 0. 

The modes in Eq. (180) are localised near n = Uc where a; = 0, in consistency 
with the linearisation in Eq. (177), and their typical width is (5a; 2± 1, or equivalently 

Sn 2± q;-^/^ . (182) 

When a <S 1 then 5n » 1 and the long wavelength approximation in Eq. (179) is 
justified. As Sn is the only scale of length in the solution (180), it must determine 
the penetration depth of the ferromagnetic magnetization into the paramagnetic 
phase, see Fig. 14. Note that this (5n ~ ^ in the general Eq. (29). 

What is more, the analytic solution (180) implies a finite energy gap in the 
spectrum of 

A = ujQ + cji = VSa ~ a^/^ , (183) 

in agreement with the general scaling predicted in Eq. (30). 

2.13.12. Inhomogeneous transition 

The general argument in Section 2.8 can be also illustrated by a solution in the 
quantum Ising chain. Consider a time-dependent transverse field gn{t) which can 
be linearised near the critical point = 1 as 

gn{t) Pi I + a{n-vt) . (184) 

A critical point at n = vt moving with a constant velocity -y > separates an 
expanding ferromagnetic phase for n < vt from a shrinking paramagnetic phase 
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Figure 15. Comparison between the quasiparticle density predicted for an inhomogeneous transition in 
Eq. (187) (solid blue), the quasiparticle density after a homogeneous transition in Eq. (121) (dotted red), 
and numerical simulations on a lattice oi N = 1000 spins (crosses) at a fixed inhomogcneity a = 2~®. 
(Figure from Ref. [45]) 



for n > vt. We want to know what is the density of quasiparticle excitations left 
behind the front in the ferromagnetic phase? 
The time-dependent Bogoliubov-de Gennes equations are 



.d , 



2gn{t)uf 



(185) 



where 
mat ion Cn 



E 



'n,m ~F "'n,m 



Unm i Vnm are Combinations of Bogoliubov modes in the transfor- 
^~d(^inm7m + <-m7m)- ^ long wavelength approximation v^^^^^n ~ 



m=Ov 



leads to approximate long wavelength equations 



idt 



[2a{n - vt)a'' + 2ifT^a„] 



(186) 



These equations were solved in Ref. [45] . Their analytic solution predicts the density 
of quasiparticle excitations 



4 ^3/4, 



rir 







c(t-q) 



when V > 2 
otherwise , 



(187) 



is the density after a uniform transition in Eq. (121) 



where n,^{TQ) = 

1/av. The analytic solution is not consistent with the long wavelength 

2. This is why the exact equations (185) were 



with TQ 

approximation in Eq. (186) near v 
also solved numerically. Figure 15 compares the analytic result (187), the numerical 
results, and the result for a homogeneous transition in Eq. (187). 

As expected from the general scaling argument in 2.8, in the Ising chain there 
is a threshold velocity u ~ 1 below which excitation of quasiparticles is suppressed 
with respect to the homogeneous KZM. Analytic solution of the Ising chain gives 
the precise value v = 2. Not incidentally, this is the fastest group velocity of 
quasiparticles at the critical point whose dispersion relation is ~ 2\k\ for \k\ <^ 
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TT. As we have seen in Section 2.13.11, the spontaneous magnetization from the 
ferromagnetic phase behind the front penetrates into the paramagnetic phase ahead 
of the front, see Fig. 14. The paramagnetic spins near the critical front are biased 
to choose the same magnetization as in the ferromagnetic phase and in this way 
excitation of kinks (quasiparticles) behind the front is suppressed. The threshold 
■0 = 2 is the fastest velocity at which the penetrating order parameter can catch up 
with the moving front. Indeed, the solution in Ref. [45] shows that when v < 2 the 
ferromagnetic magnetization penetrates into the paramagnetic phase to a depth 

/ 2 \ 1/4 

6n{v) ~ (l-x) ^^^^^ 

which shrinks to zero when u — > 2~. Above v = 2 there is no bias to suppress kink 
excitations. 



2.14. Quench across a multicritical point of the XY chain 

The XY model is a generalisation of the Ising model (92) 

H = - E + VK+i + 5<) • (189) 

n 

The Jordan- Wig ncr transformation (100) followed by a Fourier transform — 
Z^fc Cfce*'^"/-\/A^ maps the XY model to a one-dimensional quadratic fermionic 
Hamiltonian 

H = E {2 [g - {Jx + Jy) cos k] c|.Cfc + i{Jx - Jy) sin k C-kCk + h.c.j . (190) 
k 

This Hamiltonian has the general form (59) with 



Hk = 2 



g - {Jx + Jy) cos k i{Jx-Jy)smk 
—i{Jx — Jy) sin k —g + {Jx + Jy) cos k 



(191) 



The multicritical point is located at Jx = Jy = ^g, see Fig. 16, where the Hamilto- 
nian (191) is Hk = 2(7(1 — cos k)a'^ with a quasiparticle spectrum Uk = 25((1 — cos k). 
This critical spectrum is quadratic in k for small k, hence z = 2 is the dynamical 
exponent. 

We make a linear quench across this critical point along a line in the parameter 
space where Jx = ^g — ^{t), Jy = \g with the usual e{t) = t/rq, see Fig. 16. Along 
this quench line the Hamiltonian (191) is 

Hk{t) = e{t) [2{a^cosk + aysmk)] + [2g{l - cos k)a^] . (192) 

At A; = its spectrum is linear in |e|, wq = 2|e|, hence the critical exponents satisfy 
uz = 1 and, consequently, z = 2,1/ = 1/2. With these exponents the simple KZM 
estimate (17) predicts ricx ~ '^q^^^- 

Equation (192) corresponds to the general Eq. (63). The orthonormalisation that 
follows Eq. (63) leads to Eq. (65) with functions 

A(A;) = 2^(1 - cos k)\ sin k\ , (193) 
6(A;) = 25(1 -cos A;) cos A; . (194) 
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Figure 16. Phase diagram of the XY model (189) in the plane of anisotropy {J^ — Jy)/{Jx Jy) (vertical 
axis) and transverse field g/{Jx + Jy) (horizontal axis). The vertical bold lines denote Ising transitions. 
The system is also gapless along the horizontal critical line. FM^ (FMy) marks a ferromagnetic phase with 
long range order in the x (y) direction, and PM marks paramagnetic phases. The considered multicritical 
point is located at (1, 0). The red arrow indicates direction of the linear quench near the multicritical point. 

Near the Fermi point kp = we can expand A^(A;) ~ \k\^ and b'^{k) ~ /c^ and 
identify the exponents in Eqs. (71,77) as 

ZA = G , Zb = 4. (195) 

Here the inequahty za < Zh in Eq. (78) is not satisfied and the exact density 



'Q 



'I/ZA 



^1/6 



(196) 



decays with a different exponent than the 1/4 predicted by the simple KZM esti- 
mate in Eq. (17). The anomalous scaling (196) was obtained for the first time in 
Ref. [60], see also [65]. More examples of multicritical anomalous scaling can be 
found in Ref. [61]. 



2.15. Kitaev model in 2D: quench across a gapless phase 

The Kitaev model is a spin- 1/2 model on a two-dimensional honeycomb lattice in 
Fig. 17 with a Hamiltonian [64] 



H 



E 

j-|-(=even 



T „X X I 7 V 



(197) 



where j and I are row and column indices of the lattice, see Fig. 17. It is an exactly 
solvable model with a gapless phase for | Ji — J2I < J3 < Ji + ^2- The model can 
be mapped onto a non-interacting fermionic model by a suitable Jordan- Wigner 
transformation [64], 



H 



+ ^2^nO, 



n+A/2 



+ J3D, 



(198) 



where and 6.^ are Major ana fermions sitting at the top and bottom sites respec- 
tively of a bond labelled ft, vectors n = \/3ini + (^« + |j)n2 denote the midpoints 
of the vertical bonds shown in Fig. 17, and ni,n2 are integers. The vectors n form 
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Figure 17. Schematic representation of the Kitaev model. (Figure from Ref. [63]) 



a triangular reciprocal lattice whose spanning vectors are 

Ml = + \] and 

M2 = — f j- The operator commutes with H and can take values ±1 in- 
dependently for each n. The ground state corresponds to = 1 at every bond. 
Since is a constant of motion, the dynamics of the model starting from the 
ground state will never take the system out of the subspace with = 1. 

When = 1 a Fourier transform brings the fermionic Hamiltonian (198) to 
H = ^]^tp^j^H^ifjj^, where ipt = (at, &t) are Fourier transforms of a.^ and bft, and 

the sum over k extends over half the first Brillouin zone of the triangular lattice 
formed by vectors n. Here we quench 

J3(t) = e{t) = — (199) 

from t — )• — oo to t — )• +cxd. The Hamiltonian has the general form of Eq. (65) 
with 

A{k) = 2[Ji sin(feMi) - J2 sm{kM2)] , (200) 
b{k) = Ji cos{kMi) + J2 cos(fcM2) , (201) 

a{k) = 2, cr{k) = a^, and cr±{k) = . The eigenvalues of are and the 
quasiparticle spectrum = \j 4[e + hik)^ + b?(k) is gapless in the gapless phase 

where | Ji — J2I < e < \ J'A - For each e in this range there is a Fermi point kp 
such that IS.(kp) = and cj^ = 0. 

The set of all Fermi points for different e is a one dimensional Fermi line where 
sin(A;Mi) = sm{kM2)- When tq ^ 1 the excitation probability pj: in Eq. (68) is 
1 on the Fermi line and exponentially small everywhere except near this line, see 
Figure 7. Expansion ink—kp near the Fermi line gives 2:^ = 2, Zf, = 4 when Ji = J2 
and = 2, = 2 otherwise. In any case, the integration in Eq. (69) gives 



d-m 

"■ex — '^Q ^ 



(202) 
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This result was obtained for the first time in Ref. [63]. A sudden quench from a 
topologically ordered Hamiltonian to a Hamiltonian that does not support topo- 
logical order was considered in Ref. [91]. 

2.16. The random Ising chain: logarithmic dependence of excitation density 
on transition rate 

An example of the random Ising chain is [92, 93] 

N 

H = -J2is< + Jn(T-n<+l) ■ (203) 
n=l 

with periodic boundary conditions (Jn+i = <?i- Here J„'s arc random ferromag- 
netic couplings which, without loss of generality, can be assumed positive, J„ > 0, 
and 5 is a uniform transverse magnetic field. This model has two quantum critical 
points at gc = ±exp (in J„) separating a ferromagnetic phase, \g\ < g^, from two 
paramagnetic phases, \g\ > gc- The randomness is a relevant perturbation leading 
to a different universality class than the pure Ising model. 

When ensemble averaged quantities are considered, the critical exponent is = 2, 
instead of i/ = 1 in the pure case, and the dynamical parameter z diverges as 

1 



near the critical point where e = {g — gc) / gc = 0. 

In a zero order approximation [85, 94, 95], one might take first the limit e — > 
which implies z — > oo. In this limit, the general KZM estimate (16) implies ^ ~ 1 , 
i.e., neither the correlation length ^ nor the excitation density — depend 
on the transition time tq. No matter how slow the transitions is, the density of 
excitations remains the same. This approximation demonstrates that there is no 
usual power law KZ scaling in the random Ising model. However, it is too crude to 
exclude a weak logarithmic dependence. 

A more accurate result is obtained with the adiabatic-impulse-adiabatic approx- 
imation which is central to KZM [85, 94], see Section 2.3. The energy gap depends 
on e as A ~ |e|^^ ~ |e|^/l'^l while the transition rate is |e/e| = l/\t\ = l/|rQe|. 
The rate equals the gap at e when 

^ = e^/' , (205) 
TQ e 

where a 2± 1 is a non-universal parameter. When Iutq S> 1 an approximate solution 
is 

i ~ ~ (Inrgf . (206) 

This logarithmic dependence on tq is very weak as compared to any power law 
scaling. It means that no matter how slow the transition is the ensemble-averaged 
density of defects in the final ferromagnetic phase, 

nex ~ (lnTQ)-2 , (207) 

remains roughly the same. 



(204) 
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Quench time 

Figure 18. Density of Icinlcs function of transition time tq on a lattice of N = 512 sites. Circles 

are final kink densities averaged over 4 realizations of Jn- Error bars set the size of the circles. The solid 
line is the best fit nex(TQ) = g^f^^^J^ 4) where (.{tq/o) is the solution of Eq.(205). The dotted line is the 

best power lavf fit d ~ ''"q™ t° the left-most 3 data points and the dashed line is the "fit" to the right-most 
2 data points. The exponents are ui = 0.48 and w = 0.13 respectively. The exponent w = 0.48 for the 
fastest transitions is consistent with the exponent 5 in the pure Ising model. (Figure from Ref. [94]) 



The prediction (206) was tested by numerical simulation of the time-dependent 
Bogoliubov-de Gennes equations 

Here = Un,m^Vn,m are combinations of Bogoliubov modes in the transforma- 
tion Cn = X]m=o('"nm7m + ^'nm7m)- The random J„'s were uniformly distributed in 
(0, 2) implying = 2e~^. g{t) = —Qc t/TQ was quenched from a large initial Qc 
to the final 5 = where density of kinks was calculated. The initial conditions were 
stationary modes of the Bogoliubov-dc Gcnncs equations (208) at the initial large 
g corresponding to the ground state of the system. 

Results of the simulations are collected in Fig. 18. They confirm the logarithmic 
dependence in Eq. (206) for large tq. However, for relatively fast transitions that 

— 1/2 

freeze out at relatively large e, the scaling n-ex ~ Tq characteristic for the pure 
Ising model is recovered, compare Eq. (121). These transitions cannot feel the effect 
of disorder that is appreciable only when e is close enough to the critical point. 

The density of kinks scales as nex ~ (Inrg)"^. Since each kink contributes to 
the excitation energy, one might expect that the excitation energy density e at 
the final 5 = is proportional to rigx- Contrary to this simple expectation, the 
ensemble-averaged energy density was found to scale as 

e ^ (lnTQ)-(3-4±o.2) (209) 

with an exponent significantly greater than 2, see Ref. [95]. The origin of the steeper 
exponent can be traced back to the uniform distribution of J„ € (0, 2) which docs 
not exclude arbitrarily weak bond strengths J„. With increasing tq the excited 
kinks tend to localize more and more on the weakest bonds, which are the easiest 
to excite, making the excitation energy decay faster than the kink density. 

An alternative derivation of the scaling (207) from the Landau-Zcncr theory, 
taking into account distribution of minimal gaps at different values of g during 
the quench, can be found in Ref. [95]. This method originates from the theory of 
quantum annealing [12]. 
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In summary, in the random Ising model, representing the infinite disorder uni- 
versality class, density of excitations does not follow the usual KZM power law 
decay with the transition time, but the actual logarithmic dependence can still be 
obtained from the adiabatic-impulse-adiabatic approximation essential for KZM. 
This example and that in Ref. [96] suggest that disorder, when a relevant pertur- 
bation, makes excitation much easier than in the pure case, but we would need 
more evidence to support this claim as a general conclusion. 

2.17. Topological insulators: anomalous excitation of edge states 

Topological insulators are an intriguing state of matter where edge transport exists 
even in the presence of a bulk energy gap [97]. The edge states responsible for the 
transport are usually localised in the interface separating two topologically differ- 
ent insulators, the simplest example being an integer quantum Hall effect sample 
and the vacuum [98]. The topological edge states arise in a variety of systems such 
as one dimensional spin models [99] , the integer and fractional quantum Hall effect 
[100], and have been experimentally realised in topological insulators [101]. They 
also include the anomalous half-integer quantum Hall effect in the honeycomb or 
square lattice [102] and the quantum Hall spin effect [103]. Depending on the sym- 
metries of the Hamiltonian and the dimension of the system topological insulators 
can be classified in a periodic table [104]. One of the simplest systems is the one 
dimensional chain of Majorana fermions [105], where the edge states are the Ma- 
jorana fermions localised at the ends of the chain. Isolation of Majorana fermions 
would be an important step towards topological quantum computation [106]. Var- 
ious proposals include the vortex core of two dimensional p + ip superconductor 
[107], tri-junctions of superconductor-topological insulator-superconductor [108], 
and many others [109]. 
Reference [110] considers the one dimensional fermionic chain [105] 

H = ^2 (^—wojaj+i + AajOj+i — + h.c.^ (210) 

j 

where aj are spinless fermionic lattice annihilation operators satisfying {aj,a^} = 
5ij. The fermions hopping between lattice sites with the tunnelling frequency w can 
be injected or removed from the wire in the form of Cooper pairs. In the following 
we assume /j, = and imaginary A = z|A| for simplicity. 

In case of periodic boundary conditions, qn+i = ai, there are no edge states 
and the system is translationally invariant. In the quasimomentum representation 
Qj = X^fe ^it^"*^"' with k € (— 7r,7r], the Hamiltonian becomes 

with energies of Bogoliubov quasiparticle excitations = 
2^ w"^ cos^ k + I A|^ sin^ k. These are the bulk excitations of the topological 
insulator. For a fixed |A| > they are gapless at a critical point w = with the 
critical exponents z = u = 1. 
Here we consider a linear quench 



w{t) = — 



(212) 
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Figure 19. The fermionic chain (210) is mapped to a Majorana ladder by the transformation (214) which 
defines 2 Majorana fermions on each lattice site. (Figure from Ref. [110]) 



from — |A| to | A|. Equation (211) becomes a set of independent Landau-Zener prob- 
lems for each k like in Section 2.11. When tq ^ | then only quasiparticles with 
k ^ and A; vr get excited with probability pk = exp(— 7r| A^tq sin^ A:/| cos A;|). 
In the thermodynamic limit density of bulk excitations becomes an integral 



'^bulk 



^ dk 



Pk 



1 

vrlAI 



-1/2 



(213) 



when Tq S> |A|^^. This is the usual KZ scaling expected for the critical exponents 
z = 1^ = 1. 

So far everything was standard and, probably, not worth a special Subsection. 
However, as expected from a topological insulator, an open chain is different from 
the periodic chain because there are localised edge states. In case of the open chain 
it is better to introduce hermitian Majorana fermions 

%-i = ^ (e-/^at + e-/\) , c,, = -j= (e-Z^a] - e-/\) , (214) 

satisfying {ci,Cj} = 6ij, which transform the Hamiltonian (210) into 

N-l 

H = i"^[{w+ \A\)c2jC2j+l + {-W + \A\)c2j-lC2j+2] ■ (215) 

i=i 

The Majorana fermions live on a virtual ladder whose j-th rung contains two 
Majorana fermions C2j-i, C2j and corresponds to the j-th site of the original chain, 
see Fig. 19. 

At the initial w = — |A| the Hamiltonian (215) is a sum of — 1 products 
2i|A|c2j_iC2j+2, each of them having two eigenvalues ±|A|. In addition to them, 
there are two Majorana operators C2,C2n-i which, in case of periodic boundary 
conditions, would couple through the term 2z| A|c2Ar-iC2 which is missing in the 
present open chain. Their Hilbert space has two states with zero energy. The paired 
and unpaired states are shown in Fig. 20a. 

In a similar way, at the final if = |A| there are N — 1 products 2i\A\c2jC2j+i with 
eigenenergies ±|A| and two "unpaired" Majoranas ci,C2n with two zero energy 
states, see Fig. 20b. In both cases, the "paired" states are the gap-full bulk modes 
and the "unpaired" zero energy states are the edge states created by opening the 
periodic chain. 

The bulk and edge states can be written in a compact form with the help of 
the BCS state |rj) = f]^^ [uk + ^feO-fc'^l ) 1^) ' where |0) is a Fock vacuum for aj 
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(a) (b) 

Figure 20. The gap- full paired bulk states and the unpaired edge zero modes 
for the initial w = — |A| (panel a) and the final w = |A| (panel b). (Figure 
from Ref. [110]) 



and the Bogoliubov coefficients are Uk = - ^^^^,Vk = ^i^y^l + 

This state is a zero energy state of the periodic chain with a Hamiltonian H = 

Yjk>o - _7fc -) , where 7^,+ = n^Ofc + VfcO^fc 7fc - = + 

Uka}_y.. In case of the open chain and the initial w = — |A| the edge zero modes 
are C2\^) and C2N-i\^)-, and the bulk eigenstates of energies ±|A| are (±ic2j-i + 
C2j+2)\^)- At the final w = |A| the zero energy edge states are ci\Vl) and C2n\^)i 
and the bulk eigenstates of energies ±|A| are {±ic2j + C2j+i)\^)- 

We assume that the initial state at w = — |A| is one of the two unpaired edge 
states: the left zero mode |-L(— |A|)) = C2\0,) or the right one |A|)) = C2n-i\^)- 
In adiabatic transition the left and right zero modes would evolve with increasing 
w into [105, 110] 

\L{w)) oc (c2 + r C6 + • • • + C2n) \^) , (216) 
\R{w)) cc (c27v-i + r C2N-5 + ■■■ + ci) \n) , (217) 

where r = (|A| +ti;)/(|A| — w) increases with w from r = at the initial w = — |A|, 
through r = 1 at the critical = 0, to r — )• 00 at the final w = | A|. As w increases, 
the edge states get increasingly de-localised until they become uniformly spread 
along the whole ladder at the critical point. These de-localised states can be also 
written as 

\m) oc (7V+7i+)l^^) , 1^(0)) oc (7V-7i+)l^^) • (218) 

Each of these two states is an equal superposition of two gapless bulk modes from 
the positive and negative energy branch. Thus no matter how slow, the evolution 
at this point cannot be adiabatic and any further increase of w above will not 
connect these states to the edge states at the final w = |A|. An initial edge state 
will end with probability 1/2 in either the positive or negative energy band. Thus 
the probability to excite a defect is 1/2 and does not depend on tq. 
A similar analysis in Ref. [Ill] in another topological insulator - the Creutz ladder 

— 1/2 

[112] - concludes that, while the density of bulk excitations decays like Tq , the 

edge modes are excited with a probability scaling as Tq^'^^. This is much steeper 
decay than for the bulk excitations in contrast to the Majorana chain, where the 
edge modes are excited with the probability 1/2 independent of tq. Unlike in the 
Majorana chain, the edge modes of the Creutz ladder remain localised at the critical 
point. These two examples demonstrate that the edge modes, whose existence is 
ensured by topology, need special treatment in dynamical phase transitions. 
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2.18. The Lipkin-Meshkov-Glick model: KZM and infinite coordination 
number 

Most examples in this review are in one or at most two spatial dimensions, but 
in this Section we follow Ref. [113] and consider the opposite extreme of infinite 
coordination number. In this limit it is hard to identify an analogue of the KZ 
correlation length ^ but, nevertheless, the density of excitations can still be obtained 
from KZM. 

The specific model that we consider is the exactly solvable Lipkin-Meshkov-Glick 
model, 

H = - J^Y: + l^>f) -9Y.^t- (219) 

i<j 1=1 

It was introduced for the first time in the context of nuclear physics [114] and then 
thoroughly studied in a number of papers [115]. Here — t- oo is a number of spins 
in the system, 7 < 1 is the anisotropy parameter, and g is the transverse field. In 
a sense, the model is the infinite coordination number limit of the XY model or, 
when 7 = 0, the quantum Ising chain. 

The model has a second order quantum phase transition at gc = 1 with mean- 
field critical exponents. The magnetization in the x-direction (or the xy-plane in 
the isotropic case of 7 = 1) is 

m = (1-5^)^/2 (220) 

in the ferromagnetic phase when g < 1 and zero otherwise. The energy gap vanishes 
at the transition as 

A = V(ff-l)(<7-7) (221) 

for g > 1. The ferromagnetic ground state below g = 1 is doubly degenerate for 

any 7 < 1. 

A sudden quench in the Ising version of this model, when 7 = 0, was considered 
in Ref. [116]. Here we follow Ref. [113] and consider a linear quench 

g{t) = (222) 

from the ground state at g' —t- 00 to (7 = 0. In the adiabatic limit we would expect 
final magnetization m = 1 with all spins pointing in the same direction. After a 
transition with a finite rate a fraction mine of spins is reversed and the magnetiza- 
tion is incomplete 

m = 1 — mine • (223) 

This reduced magnetization costs a finite residual excitation energy per site ^jq^. 

When mine < 1 then ~ mine- 

In Ref. [113] the linear quench was evolved numerically with the results collected 
in Fig. 21a. When tq < 1 we can see saturation at mine ~ 1- These quenches 
are effectively instantaneous as there is simply not enough time to build up any 
ferromagnetic magnetization. However, when tq S> 1 then we observe a scaling 
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Figure 21. In a, residual energy per spin and mine as a function of rg in a system of N = 1024 sites. 



For slow enough quenches the data are consistent with the Tq 
as a function of tq for different system sizes ranging from A'^ 



3/2 



scaling 



In b, residual energy per spin 
32 to Af = 1024. For very large tq the 
dependence of energy on the transition rate crosses over to a steeper power law Tq^, as explained in Ref. 
[113] by incomplete LZ anti-crossings. (Figure from Ref. [113]) 



consistent with 

^ ~ mine ~ Tq'/' . (224) 

The exponent can be explained as fohows. 

Since the parity operator n^=i is a good quantum number and the initial 
ground state at large g is the even-parity state fully polarized along the z-axis, 
then we can confine to the subspace of even parity. In this subspace, when g ^ 1 
then the instantaneous gap to reverse two spins is 4(7 w On the other hand, 



m 

the gap at the critical gc = ^ scales as 



Ac ~ (225) 

when — )• oo. Consequently, the transition between the ground state and the first 
excited state can be described by an effective LZ model with a time-dependent 
Hamiltonian 



H^z = ^ / , N (226) 

and the LZ excitation probability is 

P = exp(-^A2rQ) exp(-JiV-V3,^) (227) 

Prom this probability we can read that, for a given tq, the maximal size of a defect- 
free system is A^frgg ~ "^q"^ ■ ^ similar way as in Section 2.13.2 and Ref. [36], we 
can use its inverse N^^^ as an estimate of the density of reversed spins, 

mine ^ iVf-'e ^ ^q'^' • (228) 



This estimate explains the |-scaling observed in Fig. 21a. 
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Figure 22. Schematic set up of the experiment in Refs. [2, 19], A 2D lattice is formed by overlapping two 
optical standing waves along the y-axis and the 2-axis with a Bose-Einstein condensate in a harmonic trap. 
The condensate is then confined to an array of several thousand narrow potential tubes. (Figure from Ref. 
[19]) 

2.19. Bose-Hubbard model: transition between gapped Mott insulator and 
gapless superfluid 

The non-integrable Bose-Hubbard model is one of the two paradigmatic examples of 
systems with a quantum phase transition [15], the other being the already described 
integrable quantum Ising chain. What may be even more important, the model was 
realised experimentally with ultracold atomic gases [2, 13, 19], see Fig. 22. Here 
we consider mainly its one-dimensional version 

N N 

H= - jj](4+ia, + h.c.) + - J];n,(n,- 1) , (229) 

s=l s=l 

where are bosonic annihilation operators, Ug = ala^ are number operators, and 
J is the tunnelling rate between nearest-neighbour lattice sites. For the sake of con- 
venience, we will use dimensionless units such that the on-site interaction strength 
is ?7 = 1, and assume periodic boundary conditions. In a realistic experiment with 
a few tens of lattice sites the dynamics should not be affected by the boundary 
conditions, if not to mention that the periodic boundary conditions should be di- 
rectly accessible in a ring-shaped optical lattice [117] or by painting arbitrary and 
time-dependent potentials [7], see Fig. 23. 

The model is especially interesting in the thermodynamic limit when the number 
of lattice sites N ^ oo and a number of particles is commensurate with A^, i.e., 
average number of particles per site is an integer n. In this limit it has a quantum 
Berezinski-Kosterlitz-Thouless transition [118] at Jc — Jt-"^, see Ref. [119]. When 
J > Jc the system is in the gapless superfluid phase with (algebraically decaying) 
quasi-long-range order, and when J < Jc it is in the gapfull Mott insulator phase. 
Here we consider dynamical transitions both from Mott insulator to superfluid and 
back. 

We begin with a linear ramp from Mott insulator to superfluid. 



(230) 
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Figure 23. Painting an arbitrary and time-dependent potential by a rapidly moving laser beam. In a) a 
rotating two-site potential with Bose-Einstein condensates (BEC's), and in b) a toroidal BEC which is 
adiabatically converted into 5 disconnected spots and then back into toroidal form. (Figure from Ref. [7]) 



where tq is the usual quench time. The evolution begins at t = from the ground 
state of the Hamiltonian (229) at J = which is the Mott state 

|n, n,n, . . .) (231) 

with a definite number of n particles per site. Experimentally, a perfect linear ramp 
of the tunnelling rate can be obtained by logarithmically reducing amplitude of 
optical lattice with an (optional) minor adjustment of interaction strength via the 
Feshbach resonance [55], see Figure 5, or by painting a time-dependent potential 
[7], see Figure 23. 

In a transition to the superfluid phase we are interested in correlation functions 

Cdt) = ^(V'(t)|al_,^ai + h.c.|V(t)) . (232) 

The correlations are related to the momentum distribution of atoms Uk by a Fourier 
transform Uk = exp(i/ci?)C/j. They are also good observables because Crs 

are conserved by the hopping term in the Hamiltonian (229) and by the end of 
time evolution at large J, when the interaction is just a small perturbation to the 
hopping term, they take stable final values. In particular. 



Ki = l-Ci = ^(1 - cos k)nk (233) 

k 

is a (normalized) kinetic hopping energy. Both the hopping energy in particular 
and the momentum distribution in general depend on the transition time tq. 

Since the model is not integrable, in the following Sections we review approximate 
solutions in different regimes of parameters. 

2.19.1. Slow transition from Mott insulator to superfluid 

We begin with the limit of slow transitions in infinite system. Numerical studies 
in this regime are extremely time consuming, therefore we concentrate mainly 
on KZM predictions. According to KZM, the state of the system after a slow 
transition has the characteristic length-scale ^ in Eq. (16), where z and v are 
critical exponents. For the Bose-Hubbard model the dynamical exponent z = \. 
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The Mott insulator-superfluid transition (at an integer density of particles) in a 
d-dimensional Bose-Hubbard model belongs to the universality class of the (d+l)- 
dimensional XY spin model [120]. Therefore: 

• In one dimension this mapping implies that ^ oo as in the Berezinski- 
Kosterlitz-Thouless transition and 

i - TQ. (234) 

As a result, the hopping energy of excitations should scale as 

Ki - - . (235) 

The exponent —2 means a steep dependence of the hopping energy on the quench 
time tq, which should make it easily discernible experimentally. 

• In two dimensions v = 0.67 and 

t ^0.40 ^-0.80 fnQ(!\ 

? ~ , -«-l ~ . (236) 

• In three dimensions we have the mean- field exponent v = 1/2 and 

e ~ r^/' , K, ~ Tq'/' . (237) 

This scaling was confirmed in Ref. [70] by analytic calculations in the Bogoliubov 
theory of Ref . [121]. 

The non-mean-field scalings predicted in ID and 2D have not been verified nu- 
merically, but see Ref. [122]. 

Slow transitions probe the universal critical exponents, but they require stability 
of experimental set-up over long transition time tq. This is why it may be more 
practical to experiment first with faster transitions, where the point of freeze-out e 
may be not close enough to the critical point to capture the universal scalings pre- 
dicted in this Section, but there are nevertheless (non-universal) scaling relations 
bearing witness to non-adiabaticity of the transition. In the following two Sections 
we consider two tractable regimes of parameters where such non-universal predic- 
tions can be made. In the next Section 2.19.2 we assume the density of one particle 
per site. This minimal commensurate density makes direct numerical simulations 
possible for fast quenches and limited lattice sizes. Numerical results can be fur- 
ther corroborated by approximate analytic solutions. In Section 2.19.3 we take the 
opposite extreme of large density when semiclassical methods become applicable 
and, finally, in Section 2.19.4 we analyse a reverse transition from superfluid to 
Mott insulator. 

2.19.2. Fast transition from Mott insulator to superfluid at small density 

In this Section we follow Ref. [123] and consider n = 1 particle per site. This 
minimal commensurate density allows for some exact numerical simulations. For 
fast quenches it is enough to consider short times when the wave function remains 
close to the initial Mott state (231) and can be approximated by a variational state 

b{t) (|0, 2, 1, 1 . . . ) + |2, 0, 1, 1 . . . ) + • • • ) /V2N, (238) 
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Figure 24. Short time dynamics of Ci{t). Numerics for Af = 10 lattice sites is given by the solid line 
(tq = 0.001) and the large dots (tq = 0.1). The dashed line presents Eq. (242). (Figure from Ref. [123]) 




where |ap + |6p = 1 and dynamics of a{t),b{t) is governed by 
A change of basis (a', b') = e**/^(a - b,-a - b)/\/2 yields 

Here the time-dependent Hamiltonian is precisely the Landau-Zener Hamiltonian 
in Eq. (39) but, unlike in the standard LZ model where t £ (—00,00), here the 
time evolution begins at i = from the instantaneous ground state right in the 
middle of the anti-crossing. This is the case (ii) considered in Section 2.9. 
The quantity of interest is 



An expansion of the exact solution in Eq. (56) for small tq gives 



Ci(t) _ 2 



n 3 



(242) 



to leading order in t/y/fq. The expression (242) suggests that simple rescalings 
make the dependence of rescaled C\l ^ffq on rescaled time t/ ^Jtq independent of 
the transition time tq. This prediction is confirmed in Fig. 24 where the numerical 
solutions Ci{t) for two widely different tq = 0.001, 0.1 collapse on each other. That 
a similar collapse in Fig. 25 extends also to large t comes as a bit of surprise. 
The variational Ansatz (238) is accurate under two assumptions: 

• Fluctuations of occupation numbers around the average n = 1 are small; 

• Particles can be displaced with respect to the initial Mott state (231) not more 
than to nearest-neighbour sites. 
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Figure 25. Scaling properties of Ci obtained numerically. Solid line: tq = 0.001, dots: tq = 0.03. Inset: 
the solid line is a power law fit to data for 0.001 < tq < 0.1 giving Ci(oo) = 0.501tq*^*. All data is for 
N = 10 and Jmax = 600. (Figure from Ref. [123]) 



The former assumption is self-consistent up to 



when significant number fluctuations begin to develop. This time corresponds to 



which is large, J ^ 1, for fast transitions with tq ^ 1. In a fast transition, the 
number fluctuations do not have enough time to develop before J enters the regime 
of strong tunnelling, J S> 1. In this regime, the correlators Cr are approximately 
constant, see Fig. 25, so the asymptotic Cr^oo) can be accurately approximated by 
Cji{t) at J when the number fluctuations are still sufficiently small for the ansatz 
(238) to be accurate. Indeed, Ci{t) ~ §\/7^ in Eq. (242) gives the correct scaling 
of the asymptotic Ci(oo) ~ ■s/tq, as confirmed by the numerical data in Fig. 25. 

The latter assumption, which makes all correlators Cr with R > 1 vanish, is 
relaxed in the Bogoliubov theory developed in Refs. [121] and [123]. In consistency 
with the former assumption, the theory truncates the Hilbert space to states with 
0, 1 or 2 particles per site only. The Mott state (231) is a vacuum, a site with 2 par- 
ticles is occupied by a quasiparticle and an empty site by a quasihole. The quasipar- 
ticles and quasiholes are hard-core bosons, but the hard-core constraint is relaxed 
to make their Hamiltonian quadratic. After solving corresponding time-dependent 
Bogoliubov-de Gennes equations and making the approximation Cr{oo) = CR(t) 
we obtain 



.1/2 

Q 



(243) 



J T, 



-1/2 

Q 



(244) 



Ci(oo) « 3.9 X 10-1^, 
C3(oo) « -3.5 X 10-3^^ 
C5(oo) « 1.4 X 10-^V^, 



and C21 ~ TQ for small tq <^ 1. Bogoliubov theory predicts that dominant odd 
correlations C21+1 ~ ^Jtq decay quickly with the distance R = 2l + 1. 
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2.19.3. Fast transition from Mott insulator to superfluid at large density 

In this Section we follow Ref. [34] and consider the opposite limit of large particle 
density, 1, in the initial Mott state (231). The large density regime is tractable 
by semiclassical methods, but it also makes the tight binding approximation leading 
to the Hubbard model (229) more problematic. For this reason most experiments 
simulating Bose-Hubbard model are carried out at densities close to unity [2], 
but experiments on number squeezing at large filling have been performed as well 
already at the early stage of research into possible occurrence of the Mott transition 
[13]. Apart from this, for slow enough quenches high energy degrees of freedom 
neglected in the Hubbard model are likely to remain adiabatic during the transition. 
In this sense, the model is a first low energy approximation to a ring of initially 
disconnected Bose-Einstein condensates being merged into a toroidal trap as in the 
experiment in Fig. 23 b. In this Section we will attempt to answer the question 
how does a winding number around the final toroidal condensate depend on the 
transition rate? The first experiment along these lines was done in Ref. [124] where 
3 initially disconnected condensates were suddenly merged into one by removing 
a "Mercedes" potential separating them. In many realisations of that experiment 
vortices were detected in the final condensate, see Fig. 26. 

For the sake of convenience, in this Section we rescale variables in the Bose- 
Hubbard Hamiltonian (229) so that 

N ^ N 

H = -JY^ {aUias + h.c.) + n,(n, - 1) (245) 

s=l i=l 

and the phase transition is at Jc ~ n~^. In the truncated Wigner method [125] 
employed in Ref. [34] the annihilation operators Og are represented by a complex 
field 0s, « y/n 4>s , which is normalised as ^g^i |0sP = N , and evolves with 
the time-dependent Gross-Pitaevskii equation 

= -J {(f>s+i - 2(f>s + <l>s-i) + 10.1^0. . (246) 

These approximations become accurate as n — > 00 and the critical point Jc — 
n"^ 0. Quantum expectation values are estimated by averages over stochastic 
realisations of the field 0s (i). For example, the correlation function (232) becomes 

Cr = ^{m\al+Ras + i^-c.\m) ~ ^{<P*<Ps+R + cc.) . (247) 
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Figure 26. Images of a Bose-Einstein condensate taJjen after sudden merging of 3 initially disconnected 
condensates. The dark spots are vortices created in the process (Figure from Ref. [124]). 
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Figure 27. Kinetic hopping energy Ki = 1 — Ci ~ 1 — cos A6s ft: ^ASf as a function of rescaled J/ J for 

— 2/3 

different tq is seen in A. When J <g; 1, all the plots overlap demonstrating that J = Tq is the relevant 

scale for J <g; 1. Individual plots depart from this small- J "common bundle" at J ~ 1, or equivalently 
2/3 2/3 

J/tq ~ Tq , when Ki = 1 — Ci is expected to stabilise. In B, we show Kji = 1 — Cr at J = 10 as a 
function of tq for R= 1, ...,5. Data points for each R were fitted with lines, their slopes giving exponents 
close to the ^-scaling predicted in Eq. (253) with error bars on their last digits. (Figure from Ref. [34]) 

Here the overline means average over stochastic reahsations. All realisations of 
4>s{t) evolve with the same deterministic Gross-Pitaevskii equation (246), but they 
start from different random initial conditions distributed according to Wigner dis- 
tribution of an initial quantum state. The initial Mott state (231) translates into 
initial fields 

MO) = e''- (248) 

with independent random phases Os G [— 7r,7r). The Mott state has the same fixed 
number of n particles at each site, translating into constant modulus |(/)s(0)p = 1 
and indeterminate quantum phases translating into random 6s- 

Our goal is to estimate average size of a random phase step A^s = ^s+i — ds 
between nearest-neighbour sites in a final state after transition. A random walk of 
phase around a chain of sites generates a net winding number whose average 
magnitude scales as ^/N times the average size of a phase step. We begin with a 
KZM-like estimate in a linearised Gaussian theory. 

The Gross-Pitaevskii equation (246) can be linearised in small fiuctuations S(j)s 
around a uniform large background, (/)s = 1 + 6(j)s, and dcps can be expanded in 
Bogoliubov modes as dcps = Xlfc ipkUue^^'^ + ^fc^fc^"*'^'*) with pseudomomentum k. 
For a constant J we have hk[t) = 6^(0)6"*'^'°* with Bogoliubov frequencies 

ojk = 2^J{1 - cos k) [1 + J(l - cos A;)] (249) 

and stationary Bogoliubov modes 

Uk = - Nk [1 + 2 J(l - cos h) + Wfc] , vk = Nk. (250) 

Here Mk is such that u\ — v'j, = 1. \n. the Josephson regime, when J <C 1, we 
have Vk ~ —Uk and a purely imaginary 5(j)s in = 1 + 5(j)s is a phase fiuctuation. 
However, for our random initial conditions (248), this linearisation is justified only 
for short wavelength modes of 0s, with k ~ zLtt, for whom the modes with longer 
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wavelength appear to be the (locaUy) uniform large background. From now on we 
concentrate on the short wavelength modes because they determine the variance 
of the nearest-neighbour phase step A^^. 

When k ibvr and J <C 1 then ujk ~ 2\/2J. Early in the linear quench (230) 
this uif; is so small that the early evolution of the short wavelength modes is ap- 
proximately impulse, i.e., their magnitude remains the same as in the initial Mott 
state and, consequently, A6g ~ 1 in this impulse stage. The impulse approxima- 
tion breaks down at J when the transition rate Uk/'^k equals cufe, and the evolution 
becomes adiabatic. This happens at 

J Tq^I^ (251) 

which is self-consistent with the assumption of J <C 1 when rg 3> 1. 

The crossover from impulse to adiabatic evolution at J is the key ingredient 
of KZM. In the following adiabatic evolution after J, but before J ~ 1, short 
wavelength phase fluctuations scale as b<^s ^ J"^''^ because the amplitudes do 
not change and the modes ttfe, v^, follow instantaneous stationary Bogoliubov modes 

Uk ~ —Vk ~ — 1/2(2J)^/^. Consequently, A^g has variance scaling as A0| ~ 



|(5(/>sp ~ J Given the initial condition for the adiabatic stage at J that A^f , ~ 

1, the phase fluctuations must shrink as A^| ~ A^^ _ (J/J)"^/^ ~ Tq^^^J"^/^ 
while J < 1. 

On the other hand , wh en J S> 1 then stationary modes -Ufc 1 and v^, ~ do 
not depend on J and A0| does not depend on J either. This means that A^| must 
stabilise between the regimes of J ^ 1 and J 3> 1, i.e., around J ~ 1 where it 
saturates at its final value 



Ml 



J>i 



« ^ " j~i 



^ Tq"' (252) 



scaling with a power of —1/3. This variance determines e.g. the correlator C\ in 

Kx = 1 - Ci = 1 - cosA^^ ~ Tq^^^ , (253) 

for TQ S> 1. The kinetic hopping energy per particle Ki is expected to stabilise for 
J ^ 1, when the hopping term dominates over the non-linearity in Eq. (246) and 
Ki becomes an approximate constant of motion, see Fig. 27. 

The integer winding number is a phase accumulated after N steps (divided by 
27r) 



1 ^ 

= ^T.^'^(^^+^^s) , (254) 



s=l 



where Arg(...) € (— vr, vr]. A random walk of phase, with the variance of nearest 
neighbour phase differences scaling as in Eq. (252), gives winding numbers with a 
variance 



~ AT Tq^^^ . (255) 

There are two limits where this scaling must fail. For very fast q uench es with 
TQ phases are completely random between neighbouring sites, so A^^ = 7r^/3, 
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V\/j,2 (-0.33 ±0.01) 
V\^25^ (-0.34 ±0.01) 
V\^„3 (-0.33 ±0.01) 




I 



1 10 100 1000 



1 10 100 

Figure 28. Variance of winding number measured at J = 10 as a function of tq for lattice sizes 

N = 512, 256, 128. Here point sizes equal error bars. The data points with tq > 2 and > 2 were fitted 

with the solid lines giving slopes close to the predicted — ^ in Eq. (255). is shown over a wider range 
of TQ to show the saturation for nearly instantaneous quenches, when tq < 2, and the crossover to steeper 

slope when < 2. The inset shows a rescaled W^/N for N = 512,256, 128,32,8 to demonstrate that 

Wf^ ~ Af in the KZ regime of tq > 2 and > 2. (Figure from Ref. [34]) 




and W'^ = N/12. For quenches so slow that W'^ < 1 the nature of the problem 
changes, leading to steeper falloff of W'^ with tq. Between these two limits the 
^-scaling in Eq. (255) for the winding number is confirmed by numerical results in 
Fig. 28. 

On one hand, the predicted 1/3-scaling of the winding number is limited to 
relatively slow transitions with rg ^ 1 so that J <C 1 but, on the other hand, 
the transition must be fast enough for the truncated Wigner method to remain 
accurate. This means that the crucial J in Eq. (251) must be much greater than 
the critical Jc — n~^, or equivalently tq <C n^. For density of, say, n ~ 100 atoms 
per site tq can range over 6 orders of magnitude. 

2.19.4- Fast transition from superfluid to Mott insulator at large density 

A reverse transition from the gapless superfluid to the gapped Mott insulator was 
analysed in Refs. [126, 127]. Since at the final tunnelling rate J = site occupation 
numbers are good quantum numbers, we will concentrate on time evolution 
of the occupation numbers and their conjugate phase fluctuations. An initial su- 
perfluid ground state at large J ^ n has large Poissonian number fluctuations 
~ ^/n at each site and small ~ l/v^ fluctuations of phase. By contrast, in the 
Mott ground state (231) at the final J = the number fluctuations are zero and 
the phases are random. We can expect that a non-adiabatic transition from large 
J S> n to J = will end in an excited state with finite number fluctuations. The 
conserved magnitude of these final number fluctuations will depend on transition 
rate. 

In this Section we assume large density of n ^ 1 particles per site when a polar 
decomposition Og = y^exp(i0s) with a phase operator 9s makes sense. Deep in 
the superfluid regime we can define Ug = n-\-5ns and Og = S9s, where the c- number 
n is the large particle density, while drig <C n and 66s are mutually conjugate small 
number and phase fluctuations respectively. Inserting the polar decomposition into 
an equation of motion idtUs = J(t)(as+i +as-i) + nsas, expanding to leading order 
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in the small quantum fluctuations, eliminating 69s, and making Fourier transform 
of Sus, we obtain a linearised Bogoliubov equation 



d 1 d 
di7{tjdi 



+ 4(l-cosA;)[n + J(i)(l-cosA;)]^ (5nfe = 



(256) 



with a time-dependent J{t). 

Notice that in case of constant J wc obtain the Bogoliubov spectrum of quasipar- 
ticle excitations uik = ■\/ 4J(1 — cos k) [n + J(l — cos k)\ which is identical with 
Eq. (249) up to rescaled units. Here we are interested mainly in the Josephson 
regime, J ^ Jc — ^jn and J <^ n, where the ground state is a number-squeezed 
state with a number variance 

An^g(J) ~ \/nJ (257) 
and a phase variance A^Qg ~ 1/ An^g . In this regime 

ujk « \/4n J(l - cos k) - [njf'^ . (258) 

We assume that the transition is slow enough for most Bogoliubov modes to remain 
unexcited before they enter the Josephson regime, except for a narrow range of low 
frequency modes with small which becomes negligible with increasing transition 
time. 

In the Josephson regime, a mode k becomes excited at a J when the transition 
rate J / J equals ujk- We consider two different functions J{t): 

• Like in the previous Sections on the Mott to superfiuid transition and in Ref. 
[126], we consider a linear 

J{t) = - — , (259) 

where t runs from — oo to 0. The transition rate J/J = l/|i| equals ojk in Eq. 
(258) at 

J ~ ■ (260) 

After J the evolution becomes impulse and the final variance at J = remains 
the same as in the ground state at J: 

Ar? ~ An^g(J) ~ n^'^Tg^'^ . (261) 

Self-consistency requires J in the Josephson regime or, equivalcntly, <C tq <C 
n. The transition has to be fast enough, tq <C n, for the crossover to the im- 
pulse stage to take place much above the critical point Jc where the Bogoliubov 
approximation would break down. 

• Since the tunnelling rate J depends exponentially on the strength of an optical 
lattice, we also consider an exponential 



J{t) = Jo e-*/^« , 



(262) 
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where t G (0, oo), like in Ref. [127]. The transition rate J/J = 1/tq equals ojk in 
Eq. (258) at 

J ~ n-^Q^ . (263) 

The final variance 

An^ ~ AnlsiJ) ~ (264) 

does not depend on n and its dependence on tq is steeper than after a linear 
transition, but J in the Josephson regime requires <C tq <C 1, i.e., a more 
narrow range of tq. 

The adiabatic-impulse scaling argument above has two limitations: as usual, it 
does not predict any numerical pre-factors and, what is potentially more danger- 
ous, it ignores any dependence of J on k. This is why we follow a more detailed 
calculation in Ref. [127] in the case of exponential transition. 

A new fc-dependent time variable r = — 2(1 — cosk)J{t)TQ, running form 
r — oo to r = 0, transforms Eq. (256) into a more universal 



9 , 2nrQ 



Suk = . (265) 



This equation demonstrates that, after proper rescaling of time, different modes 
Srik evolve in the same way. Their evolution changes qualitatively at r ~ 2nTQ 
when it crosses over from pure oscillations with constant frequency, ~ exp(ibir), to 
damped oscillations with increasing frequency. This crossover corresponds roughly 
to entering the Josephson regime near J 2± n. After the crossover the damping 
makes number fluctuations An shrink from the initial An 2± ^/n to its final value 
at J = 0. 
More precisely, see Ref. [127], 

6nk = e-'^"^«/'W,„^^,i/2(2zr) bk + h.c. , (266) 

where W is the Whittakcr function and bk are bosonic Bogoliubov quasiparticles. 
Here we work in the Heisenberg picture, where the state is a vacuum jO) annihilated 
by all bk- The number fluctuations at the final r = 0, corresponding to J = 0, are 

1 _ p-2iTnTQ p -, 

{{dukf) = {0\{dnkf\0) = + O te-*/-«(l-cosA:) . (267) 

Ztttq I J 

This limit is saturated and the number fluctuations freeze out when t ^ tq. Since 
the frozen {{Suk)'^) is independent of k, we obtain non-zero frozen on-site variations 

1 _ p-2nnTQ 

{<) - {us? = {{Snsf) = ^^^^ , (268) 

but vanishing two-site correlations between Sug at different sites. After a very fast 
quench with htq ^ 1 the final variance is equal to n like in the initial Poissonian 
ground state at J ^ n. These quenches are effectively impulse right from an initial 
J S> n and all the way down to J = - they are too fast for the initial state to 
adapt to the rapidly shrinking J. In the more interesting regime of utq S> 1 the 
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final frozen number fluctuations are 

An^ = (269) 
27rrQ 

in agreement with the scaling prediction in Eq. (264). The pre- factor missing in 
Eq. (264) is ^. 

When J these uncorrelated number fluctuations evolve with evolution oper- 
ator 

^-ine{ns-l)t/2 ^ ^-in'^t/2 ^-in Snst (270) 

and gradually translate into increasing phase fluctuations with a variance [127] 



A9^ = ^— . (271) 

2TTTQ ^ ' 

These phase fluctuations determine relative quantum depletion from the conden- 
sate, i.e., fraction of particles that are not condensed into a uniform condensate 
wave function with a site-independent phase. For our Bogoliubov approximation 
to be self-consistent they need to remain small, ^ 1, at least down to J when 
the ground state phase fluctuations are A6l2 ~ TQ. Thus the range of validity of 
the scaling in Eq. (269) is limited to 



n 



-1 



< TQ < 1 . (272) 



Notice that in the simple scaling argument leading to Eq. (264) the condition (272) 
is equivalent to J in the Josephson regime. 

We can conclude that the detailed calculation confirms the scaling An^ ~ '''q^ 
predicted in Eq. (264), provides the missing pre-factor of l/27r, justifies a posteriori 
ignoring the non-adiabaticity in a narrow range of long wavelength modes, and 
limits validity of the scaling to the same range of tq as the simple scaling argument 
based on the adiabatic-impulse approximation. 

The vanishing final two-site correlator between different bn^ has observable im- 
plications for a two-site correlation function Cr. In the Heisenberg picture at the 
final J = 0, the annihilation operators evolve as asit) = exp{—inst)as{0) and the 
correlation function 

CR{t) = {al^R{t)asit)) = Cr{0) {eMit{ns+R-ns)]) (273) 

decays to zero after a dephasing time ~ (An)^^. When the final J is precisely zero, 
then the decay is followed by periodic revivals every 27r, but when the final J is 
small but non-zero, then the revivals disappear after a dephasing time proportional 
to the width ~ of the quasiparticle band at small J, see Ref. [128] for more 
details. The (disappearing) periodic revivals of the initial superfluid state were 
observed in the experiment [2] after a sudden quench to the Mott phase, see Fig. 
29. 

In the following Section we continue with bosonic atoms in a quasi-lD trap, but 
this time without the tight-binding approximation leading to the Bose-Hubbard 
model (229). 
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Figure 29. Dynamical evolution of the multiple matter wave interference pattern observed after sudden 
reduction of the tunnelling rate from the superfluid to the Mott phase. Panels a...g are absorption images 
taken after hold times: Ofis, lOOfis, 150fj,s, 250^s, 350fis, AOOfis and 550fis respectively. A first distinct 
interference patten is visible (a), showing that initially the system can be described by a macroscopic 
condensate wave function. After 250^s the patten is lost (d) and after 550fis hold time it is almost 
perfectly restored. (Figure from Ref. [2]) 




Figure 30. The two scenarios considered in Ref. [42] . Here we review case a) when a one dimensional Bose 
gas is loaded into an optical lattice potential. A similar experiment was done in Rof. [8]. (Figure from Ref. 
[42]) 



2.20. Loading a ID Bose gas into an optical lattice: transition into the 
gapped phase of the sine- Gordon model 

The parameter that defines properties of a Bose gas in one dimension [129] is 
its interaction strength 7 = mg/h?p, where m is the mass of a particle, g is the 
interaction strength in the contact interaction g 6{xi — X2), and p is the linear 
density of particles. In a spatially uniform system the spectrum of excitations is 
gapless for any 7 and the low energy physics is described by the Luttinger liquid 
[130]. The parameter K of the Luttinger theory is i^T ~ 1 + 4/7 for strong and 
K ~ T^/\/^ in for weak repulsion 7 [129, 131]. In the extreme Tonks-Girardeau 
limit 7 —7- 00 of strong repulsive interactions [132], realised experimentally in Refs. 
[3-5], the particles behave like hard-core bosons which cannot occupy the same 
position in space. 

In Refs. [39, 42] the experimental set-up in Fig. 30a was considered, where by 
slowly turning on an optical lattice the initial zero temperature Bose gas is loaded 
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into a lattice potential. The lattice potential of increasing strength e{t) results in 
a sine-Gordon Hamiltonian [133] 

H = dx + {d^(t)f - 4e(t) cos(2V7r^ 0) , (274) 

where n(x) and are mutually conjugate, and K is the parameter of the Lut- 
tinger liquid before turning on the perturbation. When K < 2 the cosine term 
is a relevant perturbation opening an energy gap in the spectrum of the gapless 
Luttinger liquid. In the repulsive regime when 1 < K < 2, the massive excitations 
are solitons of mass 

A 2± €^ , (275) 

for small e, see Refs. [42] and [134]. Since in a translationally invariant system the 
solitons can be excited only as topologically trivial pairs of solitons and antisolitons 
with opposite momenta {k, —k), the spectrum of relevant excitations is 2\/ A^ + /c^. 
In the attractive regime when K < 1, the spectrum includes also breathers that 
can be interpreted as soliton-antisoliton bound states. 

Since at the critical point A = the excitation spectrum 2\/ A^ + k'^ is linear in 
I A; I, the dynamical critical exponent is 

z = 1 . (276) 

Consequently, the finite gap A translates into a finite correlation length ^ ~ A~^ 
in the ground state of the system. Given the scaling in Eq. (275), we identify the 
critical exponent 

(277) 



2-K 



in ^ - e-". 
A linear ramp of the lattice amplitude, 

e(t) = - , (278) 

with the time running from f = to t — t- oo, drives the system from the critical 
point at e = into the gapped phase with e > 0. This is an example of the general 
half-quench analysed in Section 2.6. Given the exponents z and u, we can use Eq. 
(26) to obtain the density of excited solitons and antisolitons 

nex ^ 1"^ ^ ~ Tq^/^^"-^^ . (279) 

The result (279) is corroborated by the calculation in Ref. [42] along the general 
lines of Section 2.12. 

The non-adiabaticity of condensate loading into an optical lattice was also stud- 
ied in Ref. [96] in presence of a harmonic trap potential and quasi-disorder. The 
disorder makes the transition much less adiabatic than in the pure case, in coher- 
ence with the example of the random Ising model in Refs. [85, 94, 95] and Section 
2.16. 

This Section completes our review of spinless bosons. In the next Section we 
consider spin-1 Bose-Einstein condensates. 



August 10, 2010 0:9 Advances in Physics revised2 



Advances in Physics 75 

2.21. Spin-l Bose- Einstein condensate: transition from paramagnetic to 
ferromagnetic phase 

A sudden quench in a spin-l ferromagnetic Bose-Einstein condensate was realised 
in the experiment of Ref. [135], see Fig. 31, and analysed in Refs. [136-138]. Here, 
following Ref. [138], we consider for simplicity a one-dimensional homogeneous 
condensate in a finite box. Since the following analysis is on the mean-field level, 
most conclusions can be easily generalised to higher dimensions. In dimensionless 
variables, a system placed in a magnetic field B along the z-axis has a mean-field 
energy functional 



Em = dz 

'box 



(280) 

Here = (V'l, "00) "^-i) describes the m = 0, ±1 condensate components, the wave 
function is normalised J^^^dz'^'^'^ = 1, and F^^y^z are spin-l matrices. The first 
term in (280) is the kinetic energy, the second and the fourth term describe spin- 
independent and spin-dependent atom interactions respectively, and the third term 
is a quadratic Zeeman shift originating from atom interactions with the magnetic 
field B. 

A phase diagram of the model (280) is interesting when ci < 0, like for e.g. ^''Rb 
atoms, and there is competition between the last two terms in Eq. (280). When 
we assume zero longitudinal magnetisation, = '^'^ Fz'^ = 0, whose integral is a 



36 66 96 126 156 1B6 216 ms 




Figure 31. Images of a cigar-shaped spin-l Bose-Einstein condensate after a sudden quench from the 
paramagnetic to the ferromagnetic phase. In a) transverse magnetization with magnetization density shown 
by brightness and its orientation by colour. In b) magnetization density and in c) magnetization orientation 
respectively shown in grey scale. (Figure from Ref. [135]). 
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constant of motion, then the relevant parameter is 

e = - 2 , (281) 

n\ci\ 

where n = ^^^^ is the density. At e = there is a transition between the symmetric 
polar phase for e > 0, where the ground state wave function is 

- (0,1,0) , (282) 

and the broken-symmetry ferromagnetic phase for e < 0, where the degenerate 
ground states are 



(e'^' , e^(xi+x-i)/2 2V4Ti , e'^-^ v^) • (283) 



The relative phase Xi ~X-i determines orientation of the transverse magnetisation 
{fx, fy) = {"^"^F-x^, "^"^Fy"^) in the x — y plane. The ground state is degenerate with 
respect to this orientation allowing topological defects like textures in ID, point 
vortices in 2D, and vortex lines in 3D. 

In the symmetric phase there are small thermal or at least quantum fluctuations 
around the ground state (282) which can be expanded into Bogoliubov modes [139]. 
Near the critical point the gap in the quasiparticle spectrum scales as 

A ~ . (284) 

Since by definition A ~ e'^^ for small e, we can identify vz = 1/2. On the mean 
field level considered here the critical exponent u = 1/2 and KZM predicts 

i ^ Tq^' , (285) 

compare Eq. (16), after a linear quench to the ferromagnetic phase. 

More precisely, see Refs. [136-138], what happens is that when the linear quench 
crosses the critical point e = the initial paramagnetic ground state (282) becomes 
a dynamically unstable false vacuum. The frequency of the k = Bogoliubov 
mode (contributing to the initial small fluctuations around this false vacuum) given 
by the gap in Eq. (284) becomes imaginary and the mode begins to grow quasi- 

1/3 

exponentially on the timescale of t ~ Tq in a similar way as in the classical 
transitions considered in Section 2.2. As e enters deeper into the ferromagnetic 
phase, more long wavelength modes become unstable. This linear instability grows 
until time t when its growth is halted by the quartic term in the energy functional 
(280) and the linearisation in small fluctuations around the initial paramagnetic 
state breaks down. By this time, the long wavelength modes with up to 
have become unstable and exponentially amplified with respect to the initial small 
fluctuations back in the symmetric phase. This exponential amplification promotes 
^ to the only relevant scale of length characterising the wave function ^ after the 
transition. 

In the spirit of the truncated Wigner method, in the numerical simulations of 
Ref. [138] the wave function ^ was initialised in the symmetric phase in the ground 
state (282) plus small Gaussian noise uncorrelated in space. The noise provides 
seed quantum/thermal fluctuations. Then the wave function was evolved by a lin- 
ear quench e{t) = —t/rQ to the symmetry-broken phase. Figure [139] shows two 
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(a) 



(b) 



Figure 32. Vector representation of the (magnified) transverse magnetisation 
{fx, fy) X 10^ at a) e = —0.28 and b) e = —2 during a transition with tq = 10. 
(Figure from Ref. [138]) 
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Figure 33. Snapshot of a magnetisation of the system at t = 2tq when e = —2 after a transition with 

TQ = 10. The top part shows transverse magnetisation = ^ /| + (times the box size squared L^). 

The bottom part shows the measurable longitudinal magnetisation (magnified by a factor of 10^). The 
vertical dashed lines help to see coincidences between minima of the transverse magnetisation and extrema 
of the longitudinal magnetisation. Average size of the longitudinal domains was verified in Ref. [138] to 
1/3 

scale as 5 ~ t„' . (Figure from Ref. [138]) 



snapshots of the transverse magnetisation (fx, fy) in the ferromagnetic phase. The 
magnetisation has orientation which is random, but correlated on the length scale ^. 
Topological textures in this random transverse magnetisation can be characterised 
by a winding number 

^/ dz -^Avgifx + ify) . (286) 
27r Jbox dz 

Density of textures scales as Reference [137] considers formation of topological 
textures in more than one dimnesion. 

Another quantity measured in Ref. [135] was the longitudinal magnetisation. 
Both in the experiment and in the numerical simulations of Ref. [138] the net 
longitudinal magnetisation was initially zero, J-^^^ dz fz = 0, and fluctuations of 
local magnetisation fz were small. Conservation of the net longitudinal magneti- 
sation allows for creation of a network of magnetic domains with opposite fz- The 
domains begin to form at the time ~ i of the quasi-exponential growth of the 
dynamical instability. The bottom part of Figure 33 shows one realisation of lon- 
gitudinal domains. The top part of the same figure demonstrates anti-correlation 

between transverse magnetisation niT = \ fx ~^ fy ^^'^ longitudinal magneti- 



1/3 

sation. Both are correlated in space on the length scale ^ c:^ Tq . 

The model (280) provides also another illustration of KZM in space, see Sections 
2.7 and 2.13.11. Reference [46] considers a transition in space with 



e{z) = a z 



(287) 



August 10, 2010 0:9 
78 



Advances in Physics revised2 

Dynamics of a Quantum Phase Transition and Relaxation to a Steady State 




Rescaled Position az 

Figure 34. Typical condensate magnetisation fx (order parameter) in a spatial transition (287). The solid 
line is the exact solution and the dashed line in the local approximation. In the critical regime near the 
critical point the exact solution deviates from the local approximation and penetrates into the symmetric 
polar phase to a depth of ^ ~ a^/^. (Figure from Ref. [138]) 

from the symmetry-broken phase where 2: < to the symmetric polar phase where 
2 > 0. Figure 34 shows the transverse magnetisation fx{z) which is the order 
parameter. The dashed hne is the magnetisation in the local approximation. It 
is non-zero in the broken-symmetry phase and vanishes as fx ~ ^/—^ = \foLZ 
when the critical point is approached, z — t- 0~. The solid line is the exact solution 
which deviates from the local approximation at e ~ o?!"^ and penetrates into the 
symmetric polar phase to a depth ^ ~ a~^/^. These scalings were verified in Ref. 
[46] to be accurate for a <C 1, as expected from the general argument in Section 
2.7. 

2.22. Adiabatic sweep across a gapless regime 

The main result of KZM is that, in the thermodynamic limit, a sweep across a 
second order phase transition cannot be adiabatic, no matter how slow is the 
transition rate 

S = Tq' , (288) 

because the system is gapless at the critical point and the adiabatic theorem can- 
not apply. Consequently, there is a finite density of excitations rigx (or density of 
excitation energy e) which scales with a power of the transition rate 5. The power 
generally increases with the dimensionality of the system making the density of ex- 
citations in higher dimensions decay faster in the adiabatic limit of 5 — t- 0, compare 
e.g. Eq. (17). Lower dimensional systems are generally less adiabatic because they 
have more low energy states available for excitation. As pointed out in the seminal 
paper [140], the same observations must necessarily apply to all gapless systems, 
no matter what is their dimensionality, but their non-adiabaticity can be expected 
to be more dramatic in lower dimensions. Gapless systems are quite generic as 
most systems with broken symmetries have gapless excitations like e.g. phonons, 
magnons, or spin waves. 
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In the zero temperature quantum limit, we can consider a linear sweep of a 
parameter k in a Hamiltonian, 

K{t) = Ki + dt , (289) 

starting from the ground state at an initial Ki and terminating at a final Kf. The 
final density of the positive excitation energy depends on the transition rate, e{5) > 

. According to Ref. [140], response of the system to the adiabatic sweep may fall 
into one of the following three regimes: 

A When e(6) is an analytic function of 6, then to leading order in 5 we have 

e{5) ~ (5^ . (290) 

Here the linear term is zero because the positive excitation energy cannot 
be made negative by a sweep in the opposite direction. This regime can be 
called mean field or analytic. Since e{6 — t- 0) = 0, there is adiabatic limit in 
this regime. On a more microscopic level, the origin of the quadratic scaling 
lies in the Landau-Zener transition in Section 2.9, see cases (ii,iii,iv) and 
Eq. (49) where the parameter driving the transition has a discontinuous 
time derivative just like Eq. (289). The discontinuity makes the excitation 
probability scale with a square of a transition rate Tq^ ^ 5^ instead of the 
usual exponential decay, see a more detailed discussion in Ref. [39] . 
B Even when e{d) is not analytic, it is still possible to approach the adiabatic 
limit as 

e{S) ~ (291) 

where o > but a 7^ 2. This regime can be termed non-analytic. On the 
microscopic level it originates from a combination of the S"^ scaling like in 
case A and large density of low energy excitations, see Ref. [39]. 
C Finally, in a non-adiabatic regime we have 

£{5) ~ ISl"" , (292) 

where a,b > and L is the system size. Here the adiabatic limit does not 
exist for infinite system size. Notice that it is not just the excitation energy, 
but the excitation energy density that diverges in the thermodynamic limit. 

In a system with well defined quasiparticles the non-adiabaticity of the transition 
can be alternatively classified by scaling of the density of excitations UcxiS) when 

1 (5 1 — 0. In general, the two classifications are different because the dominant low 
energy quasiparticles have a non-trivial dispersion relation. 

The existence of regimes B and C in the quantum limit was supported in Ref. 
[140] by solution of a generic low energy quadratic Hamiltonian, 

where (j)q and Hg are conjugate coordinates and momenta respectively. In the con- 
text of superfluidity compressibility. Here we choose 



Kg = K + Xq^ 



(294) 
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to cover all three regimes defined above, k is ramped as in Eq. (289). For a large 
initial Kj, the response of the system belongs to regime A with e{S) ~ S^, but when 
Ki = then 

|(5|('i+l)/4 
^ ~ Xid+l)/8 

and the response belongs to regime B. Here d is the number of dimensions. The 
density of excitations nex{S) leads to a different classification. When Ki is large, 
then the non-analytic regime B is realised for d = 1 and the analytic regime A 
when d > 2. However, when Kj = 0, then the system is non-analytic (B) when 
c? = 2, 3, but non-adiabatic (C) in ID. 

The crossover between different scalings can be further illustrated by e.g. Ref. 
[141] where a near-adiabatic parameter change within a gapless metallic or gapped 
insulating phase of the Falicov-Kimball model is studied by the dynamical mean 
field theory [142] . The excitation energy density scales with a power of the quench 
rate S whose exponent depends in general both on energy spectrum of the system 
and smoothness of the quench protocol. However, in the gapped insulator the 
exponent depends on the smoothness of the ramp only. By contrast, for a sufficiently 
smooth ramp protocol in the gapless metallic phase the exponent depends on the 
intrinsic spectrum of the system only, but this intrinsic behaviour is not observable 
when the ramp is not smooth enough. For instance, for a linear ramp there is a 
crossover to the scaling in case A. The quadratic scaling for energy density was 
also derived in Ref. [143]. 

2.23. Many-particle Landau-Zener problem: adiabatic passage across a 

Feshbach resonance 

Non-adiabatic dynamics of a single-particle quantum system can often be described 
by the Landau-Zener problem where a probability of the transition from the ini- 
tially occupied ground state to the excited state is exponentially small in the sweep- 
ing rate S. With more particles one generally encounters more LZ anti-crossings 
whose cumulative effect describes the behaviour of the system during the driving 
process. However, in many-body systems there is exponential density of energy 
levels and it is not possible to divide the evolution into a series of independent 
LZ anti-crossings: every anti-crossing takes finite time, but the frequency of anti- 
crossings increases exponentially. In experiments with ultracold quantum gases, 
where macroscopically large numbers of particles are involved, the microscopic (LZ) 
treatment is not practical, but it is often justified to use scmiclassical treatments 
like the truncated Wigner method [125]. In this context, non-adiabaticity in the 
scmiclassical models has been discussed recently in Refs. [34, 46, 136, 138, 144-152] 
and Sections 2.19.4 and 2.21. 

In this Section we consider a time-dependent Dicke model in dimensionless units 
with a time-dependent Hamiltonian 

H = -Stb^b + StS' + -^{b^S- +h.c.) , (296) 

where = Sx iSy arc spin operators, the spin 5* = N/2, 6 is a bosonic anni- 
hilation operator, and 6 is the sweep rate. When = 1 we recover the standard 
LZ model in Eq. (39) and the excitation probability P = exp(— 7r/()) is exponen- 
tially small in the adiabatic regime 5 <^1. However, here we are interested in the 
opposite extreme of a macroscopically large spin when A" S> 1. 



(295) 
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In the context of the adiabatic passage across a Feshbach resonance [55] the 
Hamiltonian (296) is equivalent to 

. N N 

H = -5tb^b+ — J] (ct ^c,t + 4^c,^) + -^5](6V,t + h.c.), (297) 

1=1 ^ i=l 

where Qo- are fermionic annihilation operators, and a ='[,■1 represents internal 
states of fermions, see the caption of Figure 5. Initially at t — >■ — oo the system 
is prepared in the ground state with 2iV fermionic atoms, and when t — t- +oo 
its instantaneous ground state becomes the state with N bosonic molecules. In 
the adiabatic limit, the time-dependent Hamiltonian (297) is meant to convert all 
atoms into molecules. However, as shown in Refs. [147, 149, 150] and outlined 
below, the sweep is never truly adiabatic, because it leaves behind a fraction of 
fermions which scales with a power of 5. 

In a semiclassical approximatio n to Eq. (296), we use the number-phase decom- 
position of the boson field, b = V Nn e**^, and the polar representation of the spin 
variables: Sz = -ycos^, Sx = ^sin^cos^, and Sy = ysin^sin,^. Combining the 
angles as ^ — </? + tt = ^, and using a conservation law Nn = 'f-(l — cos 9) to 
eliminate 9, we obtain H = — 2N (^S t n + n\/l — ncos(f)) . Finally, we rescale 
the Hamiltonian as 

H' = H/N = — 7 n — — n cos (p , (298) 

where 

7 = 2et' , (299) 

e = 5/N, and t' = Nt. Here n and (j) are mutually conjugate variables and n G [0, 1] 
is a fraction of bosons. In the following we skip the primes. 

The equations of motion that follow from the Hamiltonian (298) are n = —d^H 
and (j) = dnH. Given an initial n = n_ at t — >■ — oo, we want to find n = n+ after 
the sweep to t — ^ cx). These two asymptotic values n=p are related to an adiabatic 
invariant given by the action 

where the integral is along a closed trajectory n{t),(p{t) in the phase space. More 
precisely, / is the area (divided by 27r) of the phase space enclosed by the trajectory 
n(t), 0(t) with the convention that, when moving along the trajectory, the enclosed 
area is on the left. For a fixed 7 =Foo, the trajectories become ^ = (po—'jt, n = rizf 
and the action is 

= , 1+ = l-n+ , (301) 

compare Figs. 35a and 35c respectively, where the /'s are the shaded areas (divided 
by 2Tr). 

In the adiabatic limit of e — ;> we expect the invariant / to be conserved, = I . 
and an initial state with n_ = /_ bosons to be adiabatically converted into a final 
state with n+ = 1 — /+ = 1 — /_ = 1 — n_ bosons. However, in a sweep with a 
finite rate e we expect that / changes by A7 = /+ — /_> and the final number 



August 10, 2010 



0:9 
82 



Advances in Physics revised2 

Dynamics of a Quantum Phase Transition and Relaxation to a Steady State 




Figure 35. Phase portraits of the Hamiltoniaii (298). Panels a,b,c correspond to 7 = —6,-1.7,20 re- 
spectively. The solid dots are the fixed points, the asterisks mark the saddle points, the solid line is the 
separatrix, and the shaded areas illustrate the definition of the classical action, i.e., the area enclosed by 
the integral in Eq. (300). (Figure from Ref. [150]) 



of bosons, 

n+ = 1 - n_ - A7 , (302) 

is less than in the adiabatic hmit by A7. 

The change of the adiabatic invariant AI depends on n_. In the truncated Wigner 
method [125], the final A/(n_) in Eq. (302) has to be averaged over the initial 
Wigner function W- (n_ , 4>) . Even in the initial ground state with no bosons the 
function W-{n^) = 2N e:x.p{—2Nn^) has a finite spread n_ 2± 1/N. Thus we need 
to consider not only n_ = 0, but also small finite n_. 

The classical phase space portrait of the Hamiltonian (298) depends on a fixed 
7. When 7<— 2or7>2, there is only one fixed point where dnH = = d(pH, 
see Figs. 35a and b. At 7 = — 2 there is a bifurcation and in the range — 2 < 7 < 2 
there are two fixed points. In the same regime, there are two saddle points located 
where n = and cos^ = —7/2, see Fig. 35b. The trajectory connecting these 
two saddles, called separatrix, separates rotating from oscillating motions. In the 
adiabatic limit, the major part of the total change of the classical action AI is 
generated near the separatrix and especially near the saddle points when they 
arise during the bifurcation at 7 = — 2. 

The initial n_ is subject to quantum fluctuations ~ A'^"^ so we need to con- 
sider finite but very small n ~ N~^. Approximating 7 = —2 near the bifurcation, 
expanding to leading order in small n, introducing new canonically conjugate vari- 
ables P,Y and a new time variable s, we obtain an effective Hamiltonian 

p2 y2 y4 

H = ^-s^ + ^, (303) 

see Ref. [150] for more details. The Hamiltonian implies the Panleve equation 
^p- = sY — 2Y^ whose asymptotes were found in Ref. [153]: 

Y{s -00) = a(-s)"4 sin + In(-s) -|- const^ , (304) 

Y{s +00) = ± P(2s)~^ cos - j^^^-) ^ pQ^g^^ _ ^3Q5^ 

As s — )■ ±00 the adiabatic invariant of the Panleve equation tends to 7+ = ^ and 
/_ = ^ respectively. Using a relation between the constants a and p in Ref. [150] 
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and going back to the variables of the Hamiltonian (298) we obtain 

AI(n_) = n_ - - In f e'^""/^ - - — ln(2sin(7rC)) , (306) 

TT \ / TT 

where ^ G [0, 1) is a quasi-random variable. For a given initial n_, the final fraction 

of bosons is = 1 — n_ — A/(n_). 

There arc two physically interesting ways of taking the adiabatic limit: 

• When e — )■ before oo, or more precisely e <C we obtain the asymptote 

n+ 1 - n_ + — ln[2 sin(7rO] • (307) 

TT 

This final n+ is a random variable because the random variable n_ comes from 
the initial Wigner function W-(n-) = 2N exp(—2Nn-) and ^ from a uniform 
distribution between and 1. The last ^-term is zero on average, so it does not 
contribute to the average of n+, but it increases its variance. The amplification 
is relatively weak because e ^ 

• When N ^ oo before e — >■ 0, or more precisely <C e ^ 1, we have 




(308) 



Average of n+ over the initial distribution W-{n-) = 2N eicp{—2Nn-) is 

{n+) 1 - ^ elne . (309) 

The small initial quantum fluctuations of n_ ~ 1/N are amplified to much 
stronger fluctuations ^elne. 

In both cases we find that, unlike in the two-level LZ problem, the leading non- 
adiabatic effect scales with a power of e. There exists adiabatic limit when e — >■ 0, 
but there is no adiabatic regime below any threshold value of e. 



2.24. Summary 

In Section 2 we reviewed accumulated evidence that adiabatic dynamics in a gap- 
less isolated quantum system results in excitation which scales with a power of the 
transition rate (but with an exception of a disordered system where the excitation is 
only logarithmic in the rate). The leading motif was the adiabatic-impulse approx- 
imation essential for KZM. The main story left over some less universal features of 
the considered models, quite as well as some interesting ideas and problems that 
did not quite fit into the main narrative. Below we briefly mention three recent 
examples. 

In Ref. [154] the adiabaticity of the linear quench across a quantum phase tran- 
sition was readdressed within the adiabatic perturbation theory. There turns out 
to be an upper estimate on the excitation probability of the system in terms of an 
adiabatic dimension da- The adiabatic dimension is the dimension of the fidelity 
susceptibility of the driving Hamiltonian [155]. In the critical region, the quantum 
adiabatic dimension da = 2d + 2z — 2 Ay, where d is the dimension of the system, z 
is the dynamical exponent, and Ay is the scaling dimension of the driving Hamil- 
tonian [156]. The upper estimate implies that the excitation is negligible when the 
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linear quench time tq 2> L'^" , where L is the hnear size of the system. For instance, 
in the quantum Ising chain da = 2 and the upper bound implies adiabaticity when 
Tq ^ L?, in agreement with the exact Eq. (124) where L = N is the number 
of spins. By contrast, in the Lipkin-Mcshkov-Glick model the upper bound im- 
plies Tq ^ N^/^, while the exact adiabaticity condition in Eq. (227) is a weaker 
TQ S> A/"^/^, so in this model the upper bound is correct, but it overestimates the 
minimal tq required to make the transition adiabatic. Reference [39] introduces a 
family of generalized adiabatic susceptibilities Xm enumerated by their order m. 
Their X2 is the susceptibility considered in Ref. [154], but they argue that it is 
X4 that describes the excitation probability and density of excitations Uex for slow 
linear quenches. In the Lipkin-Meshkov-Glick model they obtain an accurate adia- 
batic condition tq 3> N"^^^ - the same as the exact Eq. (227). Both Refs. [154] and 
[39] make a very interesting connection between the adiabaticity and the fidelity 
susceptibility. 

This review is limited to isolated quantum systems, but it would not be quite 
physical to ignore the problem of adiabatic dynamics in open quantum systems. 
This problem has already been studied is some detail [85, 157-160], but it is cer- 
tainly far from being completely explored. The effect of classical and quantum 
noise acting uniformly on the quantum Ising chain was considered in Refs. [157] 
and [158] respectively. Numerical simulations for a model of local noise acting on 
the disordered Ising chain were performed in Ref. [159], and the effect of a static 
spin bath coupled locally to the ordered Ising chain was considered in Ref. [85] . The 
static spin bath was found to change the universality class of the quantum phase 
transition to that of the disordered Ising chain in Section 2.16. For weak coupling 
to the spin bath and not too slow quenches the density of excitations scales like 
in the pure Ising chain in Eq. (121), but for very slow transitions the scaling is re- 
placed by a much slower logarithmic dependence like in Eq. (206). Finally, in Ref. 
[160] the scaling theory was generalised to an open critical system and a quantum 
kinetic equation approach was formulated for adiabatic dynamics in the quantum 
critical region. It was found that for weak coupling and not too slow quenches the 
density of excitations is universal also in the presence of the external bath. Given 
the evidence at hand, we can conclude that a weak coupling to environment does 
not alter KZM for not too slow quenches, but in the adiabatic limit we have e — )• 
and the KZM happens very close to the critical point, where the system is very 
susceptible to the influence of the environment. Moreover, thermalisation dynamics 
close to a quantum critical point was considered in Ref. [161]. 

An opposite situation is considered in Ref. [162], where a central spin couples 
globally to the environment of the transverse Ising model subject to the linear 
quench like in Section 2.13. The quenched environment monitors the state of the 
central spin and its sensitivity is amplified by closeness to the phase transition. 
Decoherence of the central spin happens almost exclusively when the critical point 
of the environment is traversed and is significantly enhanced by the non-equilibrium 
dynamics. 

In the next part 3 we consider relaxation of the excited state in the last adiabatic 
stage of the evolution. 
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3. Apparent relaxation of an isolated quantum system after a sudden 
quench 

3.1. Introduction 

Once a system got excited, the question is what is going to be the fate of the 
excited state? Does it relax to any stationary state or, maybe, even thermal state? 
This question has not been thoroughly investigated in the context of adiabatic 
linear quenches, but with the exception of Section 2.13.10 above where some of 
the consequences of quantum dephasing have been worked out in the integrable 
quantum Ising chain. The reason partially is that in a linear quench it is hard 
to make a clear-cut distinction between the non-adiabatic process of exciting the 
system and the process of relaxation, because the relaxation begins already during 
the non-adiabatic excitation. 

Nevertheless, the relaxation problem itself can be formulated in a clear-cut way 
when we assume that a system is initially prepared in the ground state of an initial 
Hamiltonian Hq and then, at t = 0, a sudden quench, faster than any time scale of 
the system, is made to a final Hamiltonian H: 



In this way, the ground state of Hq all of a sudden becomes the excited initial state 
for adiabatic evolution with a time-independent H. At first sight, the discontin- 
uous sudden quench may appear very different from the smooth linear quench: a 
slow linear quench makes an effort to be adiabatic, while the sudden quench does 
not even pretend to be anything like adiabatic. In spite of this first impression, 
the difference is more quantitative than qualitative. As we know from Section 2.3, 
even in the linear quench one can often use the adiabatic-impulse-adiabatic ap- 
proximation. In this approximation the state of the system does not change during 
the impulse stage, when the parameter e in the Hamiltonian evolves from e to 
— e, but the ground state at e survives to become the excited initial state for the 
adiabatic evolution after — e. In this approximation the linear quench is effectively 
a sudden quench between e and — e. The quantitative difference is that in a slow 
linear quench the effective parameter jumps from e across a critical point e = 
to — e is small, while in a sudden quench it does not need to be small just as it 
does not need to cross any critical point. The relaxation after a sudden quench has 
been studied in many different integrable and non-integrable models. Section 3 is 
an attempt to review some well established concepts as well as more controversial 
conjectures that were formulated in this relatively new area of research. 

The first question that comes to mind is what is meant by relaxation in an isolated 
quantum system at zero temperature? After all, the evolution of the initial pure 
state is unitary so the state cannot relax to any mixed steady state. However, 
even when the state of the whole system is pure, a state of its subsystem Q, is 
in general mixed because the subsystem is entangled with the rest of the system 
Jl-*-. The state of the subsystem must be described by a reduced density matrix 
pQ = TrQx {ipl which is in general a mixed state. Thus a subsystem can be mixed 
even though the system as a whole remains pure. 

Having thus introduced mixedness, we can now attempt a definition of relaxation. 
A well defined question is: is there a mixed stationary state poo of the whole system 
such that 




(310) 



lim pa{t) = TTq± Poo 
t->oo 



(311) 



August 10, 2010 0:9 Advances in Physics revised2 



86 Dynamics of a Quantum Phase Transition and Relaxation to a Steady State 

for a finite subsystem ^11 If yes, then expectation values of local observables, with 
a finite support in the subsystem fi, relax to their expectation values in the steady 
state Poo- Even though the actual pure state of the whole system p{t) = \ip{t)){'ip{t)\ 
cannot relax, it appears to relax for local observables. There is no global relaxation, 
but there appears to be relaxation for local observers. 

The local relaxation is not the most general way a relaxation of a pure state 
can be defined. Suppose that we want the steady state poo to be a statistical 
description of a system like, e.g., a canonical or microcanonical ensemble. Even in 
classical physics a statistical description is not meant to capture all fine details of a 
state but only its most "coarse-grained" properties mainly because a too detailed 
description would not be tractable and thus not useful. From the point of view 
of observables, we expect good statistical description of a state to give accurate 
predictions of sufficiently coarse-grained observables, but we will not be surprised 
when a measurement of a complicated fine-detailed observable reveals a difference 
between the actual state and, say, a canonical ensemble. The same is true for an 
isolated quantum system: we expect that coarse-grained observables O can relax, 

lim {^{t)\Om)) = IVpoo O , (312) 

t— >-CXD 

but a measurement of a sufficiently fine-detailed observable can reveal a difference 
between the statistical poo and the actual pure state \ip{t)). In a many-body system, 
for instance, simple few-body observables may be accurately described by p^o, 
while a complicated many-body observable can reveal that the statistical poo is not 
accurate. Thus the pure state \ijj{t)) does not relax, but its approximate coarse- 
grained description appears to relax. 

Having established what is meant by relaxation in an isolated quantum system, 
we can ask what is the relaxed state Poo? A quick and in principle correct answer is: 
a diagonal ensemble. If a system does relax, then a diagonal infinite time average 
of its density matrix 

T^oo T Jo ^ 

is the first candidate for poo- Here a enumerates eigenstates of a non-degenerate 
H and 

V. = \{a\m)f = Kc^mm' m) 

is a conserved probability that the isolated system is in the eigenstate a. Indeed, if 
a system relaxes to a steady state, then expectation values in the steady state must 
be the same as in the infinite time average p. Unfortunately, the diagonal ensemble 
p may be disappointing as a statistical description, because it contains a lot of 
information about the initial state \ip{0)) encoded in the microscopic initial proba- 
bilities Pa- Thus p itself may be not the desired tractable statistical description poo, 
but it is a good starting point to look for a more tractable poo as a coarse-grained 
description of p. 

By analogy to classical statistical physics, one can hypothesise that the steady 
state Poo is the state of maximal entropy subject to the constraints imposed by 
integrals of motion 7^, where the integrals Im commute with H and between 
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themselves. The entropy is maximised by a generaUsed Gibbs ensemble (GGE) 

pGGE ^ exp^-^A^I^^ , (315) 

where the numbers A^^ are fixed by the conserved expectation values Trp^^^I^ = 
(V'(0)|/m|^(0)). The special case of only one integral h = H is the canonical 
ensemble. However, unlike a classical system, any quantum system has as many 
integrals of motion as the dimension of its Hilbert space. Indeed, we can always 
choose the integrals to be projectors = on the eigenstates of H, or as 

integer powers of a finite Hamiltonian = H'^ with m = 1, dim(i7). With 
the former choice the GGE becomes the microscopic diagonal ensemble p, which 
may be too accurate to be tractable. The latter choice is also equivalent to the 
diagonal ensemble p, but it may be a better starting point for further coarse- 
graining approximations, at least in the model considered in Section 3.9 below, 
where the sum over i?™ can be accurately truncated to a few lowest powers. Thus, 
unlike in a classical system, the problem with GGE is not how to find any integrals 
of motion, but how to find a subset of the abundant integrals of motion leading to 
a sufficiently accurate but still tractable GGE. Nevertheless, there is a wide class 
of quadratic bosonic/fermionic Hamiltonians, where the number operators of 
bosonic/fermionic quasiparticles provide such a small subset, see Sections 3.3 and 
3.4 below, where these non-interacting systems are shown to relax to GGE for local 
observables O. The absence of relaxation to a thermal state was experimentally 
observed in the seminal experiments on one-dimensional hard-core bosons [3, 4], 
see Fig. 3.1. 

Thus the integrable quadratic Hamiltonians relax locally to a GGE, but what 
happens with an interacting system that cannot be mapped to any non-interacting 
Hamiltonian? It turns out that in certain circumstances a non-integrable system 
may relax to a microcanonical ensemble. This scenario requires two conditions to 
be met: the eigenstate thermalisation hypothesis (ETH) proposed in Ref. [163] and 
a narrow energy distribution pa of the initial state in the eigenbasis of the final 
H. The narrow distribution often happens to be the case, see the examples in 
Sections 3.9 and 3.10. The ETH asserts that any few-body observable O has the 
same expectation value in each eigenstate \a) of a many-body Hamiltonian H in 
a narrow energy window AE. This not-quite-intuitive hypothesis is shown to be 
the case in the examples of non-integrable systems considered in Section 3.10 while 
it is not true for the integrable models considered there. When we assume these 
two conditions, then for any few-body observable the diagonal ensemble (313) is 
indistinguishable from a microcanonical ensemble 

Poo = M (31^) 

a, \E„-{H)\<AE 

where AE is a narrow energy window around average energy {H). A coarse-grained 
few-body observable has the same expectation value in p as in because, thanks 
to ETH, it actually has the same expectation value in any eigenstate in the energy 
window. In fact, assuming that ETH is exactly true, we can push it to its logical 
limit and claim that for any few-body observable the state 

Poo = , (317) 

with la) being any individual eigenstate in the energy window AE, is as good an 
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approximation to the actual diagonal ensemble p as the microcanonical ensemble. 
Indeed, since all expectation values {a\0\a) in the window are the same, there 
is no need to average over the window as in the microcanonical ensemble (316). 
Thus each individual many-body eigenstate in the window provides a thermal (mi- 
crocanonical) average for a few-body observable O. This thermal character of the 
individual eigenstates is hidden by their initial coherent superposition \t/j{0)), but 
it is revealed after the initial phase coherence is destroyed by relaxation or, more 
precisely, it appears to be destroyed for the simple (coarse-grained) few-body ob- 
servables. 

The ETH is argued to hold in non-intcgrablc quantum systems, while the GGE 
often is an accurate description of the relaxed state in integrable quadratic systems, 
see Sections 3.3, 3.4, 3.5, 3.6, and 3.7 below. This is similar to classical physics, 
where integrable systems constrained by their conservation laws do not thermalise, 
but non-integrable chaotic systems do thermalise [164]. However, this mechanism 
of quantum thermalisation is qualitatively different than that of classical ther- 
malisation. According to ETH, each many-body eigenstate is a thermal ensemble, 
while the ergodicity means that classical thermalisation requires probing essentially 
all states on a manifold of definite energy. A classical integrable system becomes 
chaotic when a non-integrable perturbation is stronger than certain threshold. A 
good example is the classic Fermi-Pasta-Ulam numerical experiment, where the 
perturbation is too weak to make the system chaotic [164]. An interesting ques- 
tion, addressed in Section 3.10, is if a similar threshold exists in a quantum system. 
Numerical simulations in small systems suggest that turning-off a non-integrable 
perturbation to an integrable model results in a smooth departure from ETH, but 
it is not clear what happens in the thermodynamic limit. 

The microcanonical steady state required both ETH and a narrow initial energy 
distribution Pa- Since in an isolated system the distribution is constant, there is no 
way it could relax to a canonical distribution pa ~ exp {—f3Ea) unless it happened 
to be canonical from the very beginning. Thus there is no global thermalisation, 
but one can still try to identify local or coarse-grained observables O for which the 
global state appears to relax to the canonical ensemble. For instance, in a system 
of many identical particles one is often interested in a single particle momentum 
distribution ra^. This interest is well motivated by current ultracold atoms exper- 
iments, where it is routine to measure the momentum distribution in an atomic 
cloud expanding after opening a trap, see the example in the right panel of Fig. 
3.1. Thermalisation of is inhibited by integrability of the system like e.g. for 
quadratic Hamiltonians which relax locally to a non-thermal GGE. It can also be 
inhibited for particles interacting by two-body interactions in one dimension, where 
the conservation of momentum and kinetic energy in a collision of two particles 
does not allow for any changes in the momentum distribution. A spectacular recent 
quantum Newton cradle experiment [4] demonstrates this inhibited relaxation in 
a quasi-one-dimensional regime, see Fig. 3.1. Integrability is also responsible for 
such phenomena as the ballistic character of transport [165] or the meta-stability 
of solitons [166]. 

Another interesting question, in addition to the nature of the final relaxed state, 
is the dynamics of relaxation to this state. Since the ground state of an initial Hq 
may contain highly excited eigenstates of a final Hamiltonian H, the answer to this 
question is not universal, because the usual universal low energy theories may not 
describe the relaxation accurately. This non-universality is illustrated by examples 
in Section 3.11, where it is shown that the relaxation can be actually the fastest 
when H is either at a critical point or close to it. This is contrary to the usual 
notion of critical slowing down derived from a universal low energy theory. In spite 
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Figure 36. Left panel: The Newton cradle experiment in Ref. [4] . An atomic 
cloud of effectively one-dimensional bosons confined in a harmonic trap is ini- 
tially split into two clouds with opposite momenta. Then the two halves keep 
oscillating in the trap and pass through each other without thermalisation of 
momentum distribution, see the right panel. (Figure from Ref. [4]) 
Right panel: The red curves in panels a,b,c are the actual momentum distribu- 
tions measured for different interaction strengths of one-dimensional hard-core 
bosons in Ref. [4]. The dashed line in panel c is the best Gaussian fit to the 
actual distribution. To the extent the actual distribution docs not conform to 
a Gaussian, the atoms have not thermalised. (Figure from Ref. [4]) 



of that, there are some remarkable universal features, like the quasiparticle horizon 
effect described in Section 3.2. 

Section 3 is organised as follows. In Section 3.2 we review the quasiparticle hori- 
zon effect. In Section 3.3 the generalised Gibbs ensemble is introduced for quadratic 
Hamiltonians, and in Section 3.4 it is shown how a quadratic Hamiltonian prepared 
in an initial Gaussian pure state relaxes locally to GGE. The following Sections 3.5, 
3.6, and 3.7 describe examples of quadratic systems relaxing to GGE, and discuss 
applicability of the GGE to different observables. In Section 3.8 we consider the 
non-integrable Bose-Hubbard model in different limits when it is close to integra- 
bility and GGE is an accurate steady state, but we also cite some evidence that it 
thermalises away from these integrable limits. Section 3.9 provides an example of 
a system integrable by Bethe ansatz, where an accurate non-thermal GGE can be 
constructed out of integer powers of the Hamiltonian. In Section 3.10 we review the 
microcanical ensemble in connection with the eigenstate thermalisation hypothe- 
sis (ETH). Finally, in Section 3.11 we mention some aspects of the dynamics of 
relaxation. We conclude in Section 3.12. 
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3.2. Quasiparticle light cone effect in dephasing after a sudden quench 

At t = a system is prepared in the ground state j-i/'o) of a Hamiltonian Hq and 
then, after a sudden quench, it evolves unitarily with a different Hamiltonian H 
at t > 0. The initial Hq is non-critical and it has a finite correlation length 
and a finite gap Aq. The sudden quench is faster than A^^. The question is how 
do the correlations evolve in time? As shown in Ref. [167], the answer to this 
question is quite universal when If is at a critical point. Particularly powerful 
analytic results can be obtained in one dimension when the dynamical exponent is 
2: = 1 or, equivalently, there is a linear quasiparticle dispersion relation = v\k\ 
with quasiparticle velocity v. This 1 + 1 dimensional problem can be described 
asymptotically by a boundary conformal field theory [168], where the initial state 
is a boundary condition. 

The results from the conformal field theory suggest a simple physical picture. 
The initial state |V'o) is a highly excited state as compared to the ground state of 
the final Hamiltonian H. It is a source of quasiparticle excitations. Quasiparticles 
originating from closely separated points, within the correlation length ^0 in the 
ground state of Hq, are quantum entangled. Once they are emitted, they behave 
semi-classically travelling at speed v. This simple picture has several important 
implications: 

• Incoherent quasiparticles arriving at a given point from well separated sources 
cause relaxation of most local observables at this point to their expectation values 
in the ground state of H. There are exceptions like the local energy density which 
is conserved. This relaxation is exponential ~ cxp (— 7r.Ti;t/2ro), where tq ~ Aq ^ 
is not universal, but x is the bulk scaling dimension of a given observable. 

• Entangled quasiparticles arriving at the same time t at points with separation 
R ^ ^0 induce correlations between local observables at these points. Since 
they travel at a definite speed v, there is a sharp light-cone effect, i.e., the 
connected correlation functions do not change significantly from their initial 
values until time t ~ R/2v. The light-cone effect is rounded off over the region 
t — R/2v ~ Tq because quasiparticles remain entangled over this time scale. 
After t ~ R/2v the connected correlations rapidly relax to time-independent 
values. At asymptotically large separations R, but well within the light cone 
where R <^ 2vt, they decay exponentially cxp {—TrxR/2vTo). This decay is 
qualitatively different from the power law decay in the ground state of H which 
is at the critical point. 

• Entangled quasiparticles arriving at the same time t 2± R/2v at points separated 
by i? ^ ^0 induce entanglement between these points. The entropy of a sub- 
system of size L, which is dominated by pairs of quasiparticles entangled across 
the boundary of the subsystem, initially grows linearly with time as more and 
more quasiparticles cross the subsystem boundary until it saturates at t ~ L/v. 
In the thermodynamic limit, in a system with periodic boundary conditions the 
saturation time is t ~ L/2v, and in an open system the entropy saturates at 
t L/v. In a conformal field theory the saturation time is rounded off on the 
time scale Aq 

• In the thermodynamic limit, the saturated entropy of the subsystem of size L is 
linear in L because the number of quasiparticles in the subsystem is linear in L 
and, after the saturation at t = 0{L/v), every quasiparticle which happens to 
be in the subsystem is more likely to be entangled with a quasiparticle outside 
than a quasiparticle inside the subsystem. Thus, there is no area law for the 
entropy, so the relaxation process cannot be efficiently simulated by, say, the 
DMRG algorithm [169] on a classical computer. 
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The light cone effect is even more rounded off when the critical H is a lattice 
system. The saturation begins at to ~ R/2vq, where -uq = ^(/c = 0) is a long 
wavelength quasiparticle group velocity, but after the entropy (or correlation 
function) keeps varying slowly with time because on a lattice there are quasiparti- 
cles slower than vq. 

Beyond the conformal field theory, the light cone effect was confirmed in the 
quantum Ising model [167], the chain of coupled harmonic oscillators [167], the 
Heisenberg chain [170], and a number of other systems, see Ref. [171] and the 
following Sections. It was also confirmed away from criticality, where vq becomes the 
maximal group velocity of quasiparticles, see e.g. Ref. [172]. In general, a maximal 
velocity of propagation of information in non-relativistic systems is guaranteed by 
the Lieb- Robinson theorem [173]. 

By contrast, in a disordered Heisenberg chain the entropy seems to grow loga- 
rithmically with time, at least as far as it can be inferred from numerical results 
[170]. This effect cannot be explained by diffusion of quasiparticles replacing the 
ballistic motion in the pure case. 

The light cone effect is the dominant feature of quantum relaxation after a sudden 
quench. The question what is the final relaxed state is the topic of the following 
Sections. 



3.3. Generalised Gibbs ensemble (GGE) for a quadratic Hamiltonian 

Integrals of motion la mutually commute, [la,!^] = 0, and commute with the 
Hamiltonian, [Ia,H\ = 0. In a unitary evolution of an isolated quantum system. 



\m) 



-itH 



IV'(O)), all moments of the integrals of motion are conserved, 

{m\ii'ii'---\m) = wo)iir/r...iV'(o)), 



(318) 



because each product /[^/2^ ... is also an integral of motion. 

When degeneracies in the spectrum of the Hamiltonian can be ignored, then 
the infinite time average of the density matrix of the system p{t) = \ip{t)){'4){t)\ is 
diagonal in the eigenbasis of the integrals of motion. 



1 

P = lim Tf; / dt p{t) 



X] Pni,n2,...\'rii,n2,...){ni,n2,...\ , (319) 



ni,n2, 



where Ia\'ni,n2, ■ ■ ■) = na\'ni,n2, ■ ■ ■)■ The diagonal p conserves all moments 



'1 -'2 



Trpll^i;^ 



1 r 

lim - dtTTp{t)Il^i;^... 

lim i rdt {m\ii^i2'---\m) 

^oo I Jo 



= {m\i['ii'---\m) , 

as expected from a candidate for a stationary state of an integrable system. 
The diagonal average p can be rewritten in an equivalent form 



(320) 



p = M exp 



a a</9 ot<P<n 



(321) 
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where the A's and the normahsation are chosen so that Tr p I]^l2^ ■ ■ ■ = 
(V'(0)|/p/2' ••• 1^(0)) and Tip = 1. The density matrix (321) has a form of a 
most general Gibbs ensemble and as such it is manifestly the density matrix that 
maximises the von Neumann entropy subject to the constraints imposed by all the 
conserved moments (V'(0)|I[^/2^ . . . |V'(0))- However, this exact representation of p 
may be not a useful coarse-grained statistical description because it contains as 
much information as the diagonal ensemble (319) itself. 

Nevertheless, the exact form (321) is a good starting point for more tractable 
approximations. The crudest but still non-trivial approximation is known as a 
generalised Gibbs ensemble (GGE) [174] 

p ^ M exp 

where the A's are fixed by conservation of the first moments only: 
Tr p la = {tp{0)\la\ip{0)). The mixed state pgge is the maximal entropy state 
subject to these constraints. The drastic approximation has a price of course: there 
are observables O for which Tr pcGE O is a very bad approximation to the exact 
TrpO. However, this is expected in statistical physics, where a sufficiently complex 
observable O can reveal a difference between a coarse grained state like pcGE and a 
microscopic state like p, but there is practically no difference for sufficiently coarse 
grained observables. In the following we attempt to figure out what are the condi- 
tions for an observable to be "sufficiently coarse grained" in case of the generalised 
Gibbs ensemble. 

To be more specific, we assume that the integrals are numbers of bosonic or 
fermionic quasiparticles [175, 176], /q = na = "ytja- Since /q,'s mutually commute 
we can factorise (322) 



PGGE = 




(323) 



This diagonal matrix is a classical ensemble for the set of independent random 
variables with a product joint probability distribution 

Pn^L- = llJ^aewi-Kna] ■ (324) 

a 

In the GGE the occupation numbers Ua are independent and, in particular, un- 
correlated. For an observable O to have an accurate expectation value in GGE the 
observable cannot depend on correlations between different n^'s. 

In case of fermions, when Ua is either or 1, this is also a sufficient condition 
because the expectation value {ua) plus normalisation determine completely the 
probability distribution for Ua- The distribution can be written in e.g. the GGE 
form Ma exp [— Aana] with the coefficients Ma and Aq determined by the two con- 
straints. 

In case of bosons, when = 0, 1, 2, . . . , the expectation value (n^) and normal- 
isation are not sufficient to determine the probability distribution for n^. Thus the 
GGE ansatz J\fa exp [— A„n„] may have wrong variance and higher moments even 
though it has the correct expectation value. 

Thus we can conclude, that GGE is an accurate statistical description for ob- 
servables which do not depend on correlations between quasiparticle occupation 



= PGGE , 



(322) 
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numbers and, in case of bosons, higher moments of individual occupation num- 
bers. This is the meaning of the mathematical conditions given in Ref. [175]. 

The GGE is an accurate statistical description of p for observables O that do 
not depend on correlations between tt-q,. Some correlations between n^'s are al- 
lowed provided that their contribution is negligible in the thermodynamic limit. A 
good example are ID hard core bosons which can be represented equivalently by 
non-interacting fermions. When perturbed by an alternating potential, as in Ref. 
[174], the initial state contains correlations between pairs of occupation numbers 
Uk and Uk+TT only. These sparse correlations make negligible contribution to lo- 
cal observables restricted to a finite region of real space. Another example is the 
Luttinger model equivalent to non-interacting bosons, see Ref. [177] and Section 
3.6. A quench of interaction strength in this translationally invariant Hamiltonian 
prepares an initial state with correlations between pairs of bosonic occupation num- 
bers n/j and n_fc only. Again, in the thermodynamic limit these sparse correlations 
have negligible contribution to local observables. 

Moreover, with an additional assumption that the initial pure state is Gaussian 
and the Hamiltonian is quadratic in bosonic/fermionic annihilation operators, GGE 
can be shown to be the unique steady state for local observables. This and the 
dynamics of local relaxation to this stationary state is the subject of the next 
Section. 



3.4. Local relaxation to GGE after a quench in a quadratic Hamiltonian 

In the last Section we discussed for what observables the time-averaged diagonal 
density matrix p in Eq. (319) can be accurately represented by the generalised 
Gibbs ensemble (GGE) in Eq. (322). Our interest in the time-averaged density 
matrix was motivated by the fact that if the long time limit of an expectation 
value Tr p(t) O exists, then it is equal to the expectation value in the time-averaged 
density matrix, 

lim Tr p{t) O = Tr p O . (325) 

t— >-oo 

However, in an isolated quantum system a limit p(t — )• oo) does not exist ^ , so the 
existence of the limit on the left hand side of (325) should not be taken for granted 
for an arbitrary observable O. This is why in this Section we investigate a more 
refined question if we can define a subset of observables O for which 

Um Tr p{t) O = Ti p^^^ O ? (326) 

t— >oo 

In other words, we are looking for a subset of observables for which a pure state 
p{t) appears to relax to a mixed steady state p*^^^. 
Following Ref. [176], we consider a general quadratic lattice Hamiltonian 



(327) 



where are bosonic or fermionic annihilation operators. The Hamiltonian is di- 



^ Except for the trivial case of a time- independent p diagonal in the eigenbasis of the Hamiltonian. 
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agonalised to if = Ylk ^klklk = Ylik ^klk by a Bogoliubov transformation 

^ = XI {^^klk + v*^,klk) ' (328) 
fc 

where the index k enumerates quasiparticle states. The quadratic Hamiltonian 
(327) has a Gaussian BogoUubov vacuum ground state annihilated by all quasi- 
particle annihilation operators Any Gaussian state p is fully determined by its 
quadratic correlators 

Q'mn ~ Tr p Cm^^n i Pmn ~ Ti' p CmCn . (329) 

Here we assume that the initial pure state p(0) is Gaussian. It can be prepared as 
a ground state of a different quadratic Hamiltonian Hq before the Hamiltonian is 
suddenly quenched to the final H. The quadratic H evolves p(0) into a pure state 
p(t) which is also Gaussian. 

The lattice can be divided into a finite subsystem Q and its environment ri-*-. We 
are interested in local observables O with a finite support in J7. To find expectation 
values of the local observables it is enough to know a reduced density matrix 

Pnit) = Ttn^ pit) . (330) 

We want to show that, under certain conditions, the reduced density matrix tends 
to a steady state, 

lim pn{t) = Ti^. p^GE , (331) 

t— >-oo 



obtained by reduction of a Gaussian GGE 



^GGE 



= ATexp 



-^^'klhk 



(332) 



Here A's are fixed by initial conditions Trp'^'^^7]^7fc = Trp(0)7j!7fc = for the 

integrals of motion 1^ = 7^7*: • If the limit exists, then for all local observables O 

in the finite the state appears to relax to GGE. This is apparent local relaxation 
which should not be mistaken with any global relaxation of the pure state p(t). 

When p is a Gaussian state fully characterized by the correlators (329), then the 
reduced density matrix p^ is also a Gaussian state fully characterised by a subset 
of the same correlators (329) with indices m, n restricted to the subsystem fi. Thus 
all that we need to show for a Gaussian initial state is that 

lim amn{t) = a^™ , lim /3mn{t) = , foT m,n en. (333) 

t-^oo t—^oo 

Here 

amn{t) =^Tr p{t) (umklk + V*^_kl'^_f}j (<fc,7^, + Vn-k'l-k') , (334) 

Pmnit) =Y^Tl p{t) (umklk + V*^-kl-k) {^nk'^k' + <-fc'7lfc') , (335) 



August 10, 2010 0:9 Advances in Physics revised2 



Advances in Physics 95 

while the correlators in the GGE are "diagonal" in k 

(^mn^ = [«mfc<fc(l ± rik) + V*^,,Vnknk] , (336) 
k 

Pmn^ = Y [UmkVnki^ ± njfc) + V^f,Unknk] , (337) 
k 

with the upper/lower sign for bosons/fermions. They follow from diagonal expec- 
tation values in GGE: 

Trp^^^'jklk' = , (338) 
IVp^^'^TikTi' = h,k'il ±nk) . (339) 

Equations (333) are satisfied if a contribution of quasiparticle correlators 

IV p{t) 7fe7fe' = e-^*^-'' TV p(0) jklk' , (340) 
Tr p{t) 7ik7iE' = e-'*^^"-'"''^ TV p(0) 7jk7j, for k ^ k' (341) 

to correlators (334,335) vanishes when t — t- oo. 

This seems to be quite generic. Indeed, given that the indices m, n are bounded to 
a finite J7, we need to consider a finite number of Bogoliubov coefficients Ufnk,Vmk 
enumerated by m € $7. In the thermodynamic limit of infinite lattice, when t — t- oo 
then the phase factor e^'^K^k+i^k') (340) oscillates with k and k' fast enough to 
average out to zero the contribution of this term to the infinite sum over k,k' . The 
same is generically true for Eq. (341) when k ^ k' . More rigorous mathematical 
conditions, and some important exceptions, can be found in Ref. [176]. 

After the dephasing in a finite subsystem O is completed, then for local observ- 
ables in J7 the quadratic quasiparticle correlators appear to be equal to quadratic 
correlators (339) in GGE (332). Since GGE is Gaussian and quadratic correlators 
determine a Gaussian state uniquely, then GGE is the unique (apparent) stationary 
state for local observables in O. 

Returning to the issue of relaxation, we can refine the above argument in a trans- 
lationally invariant case, when Vmn = Vm-n and Wmn = W„i-n in the quadratic 
Hamiltonian (327) and the Bogoliubov modes in Eq. (328) become Umk = e*^"*^^ 
and Vjnk = e^^^'Vk with definite quasimomenta k (we suppress vector notation 
in more than one dimension). A quench from a different translationally invariant 
quadratic Hamiltonian prepares an initial state p(0) such that 

Tr p(0) ^klk' = S.k,k' Afe , Tr p(0) 7^7^, = 6k,k' (1 ± n^) • (342) 
After simple algebra using the symmetry Uk = oj-k we obtain 

—e'^^^-^) (e-2^*-'=Ufe«_feAfe + e2^*-'=^;M,A^) (344) 
In a small neighbourhood (k^ — 5k, ko + 6k) of a generic ko we can linearise in 
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{k - ko) as Uk ~ ojko + Vka{k - ko), where Vk^ = ^(^o) is a group velocity ^ . 
Contributions to the above integrals from this neighbourhood are proportional to 



/Sk 
-Sk 



d{k ko) ^i(k-ko)\m-n±2vk„t] ^ ^3^5^ 



provided that 5k is much less than the shortest scale on which the functions 
Uk,Vk,nk can vary in k. If U}^^vi.,n]. are not singular ^ , then for t — >■ 00 this 
integral decays like 0{Sk/t) except when 

|m — n| «i 2 Vkgt . (346) 

The distance 2vkot is a separation of a pair of quasiparticles with opposite momenta 
(A:o, —ko) created by the quench, see the quasiparticle horizon effect in Section 3.2. 
When t is long enough this separation becomes much longer than the subsystem $7 
and the exception (346) does not apply to any m,n E ^l. This is another example 
of the quasiparticle horizon effect in Section 3.2. 

Thus the relaxation in Eqs. (333) seems to be generic and the quadratic corre- 
lators tend to the quadratic correlators in the GGE (332). Since both p^it 00) 
and Trf^-L/?*^^^ are Gaussian, the equality of their quadratic correlators implies 
that both states arc the same and, in particular, all local observables O G Q (and 
not only quadratic correlators) have the same expectation values in the two states. 
The GGE (332) is determined uniquely by the quadratic correlators plus the re- 
quirement that GGE is Gaussian. 

The Gaussian GGE (332) is determined by the initial quasiparticle occupation 
numbers n^, so it is not surprising that it correctly "predicts" the conserved n^. The 
free lunch is that after the local relaxation GGE also correctly describes all local 
observables O e Q. Thus we need to invest all to predict all local observables. 
This is lesser return than in the Gaussian canonical ensemble 



^canonical ^ J\fexp[-l3H] = M exp 



(347) 



which is a special case of GGE (332) with = Pujk and only one adjustable 
parameter /3 that we need to fix. However, this is the price for the integrability 
of the quadratic Hamiltonian. In a sense, in the GGE (332) each non-interacting 
quasiparticle has its own inverse temperature f3k = \k/oJk- 

We can conclude that in general a quadratic Hamiltonian initially prepared in 
a pure Gaussian state relaxes locally to GGE. More precisely, expectation values 
of observables within a quasiparticle horizon relax to their expectation values in 
GGE. 



3.5. GGE and the transverse quantum Ising chain 

The argument in Section 3.4 applies to the transverse field quantum Ising chain 
(92). A Jordan- Wigner transformation maps the Hamiltonian (92) to a transla- 
tionally invariant quadratic Hamiltonian (104). A sudden quench from an initial 



^Whcn Dfcj-i = 0, then wc can expand = ej-^ + A{k — feo)^ and, when Uh,Vk,nk are not singular, the 
integrals decay like 0{5k/\/i) for t — > 00. 

-"-The Mfc,ff;,nfc cannot be singular for fcrmions, when they are bounded by the constraints |itfep + = 
1 and < rife < 1, but they can in principle be singular for bosons because the bosonic constraints 
— kfcP = 1 < riA: < 00 do not provide any upper bounds. 
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transverse field gi to a final gf > prepares an initial Gaussian state which later 
relaxes locally to the GGE (332), where is an annihilation operator for a Bo- 
goliubov quasiparticle at the final gf, see Eqs. (108,109,110,111), and are such 
that Tr pGGE ^t^^ ^ p(o) ^l-fk. 

When expressed in terms of fermions, the quantum Ising model relaxes locally to 
GGE. However, since the Jordan-Wigner transformation (100) between the original 
spin operators and the fermions is non-local, there are local spin operators which 
are not local in fermionic representation, the simplest example being o"^. Thus in 
general only spin observables with a finite fermionic support relax to GGE. 

The dephasing to GGE (332) was studied in Refs. [83, 178], although without 
any explicit reference to GGE for good chronological reasons. In particular, in Ref. 
[83] ferromagnetic correlation functions C|j^ = Tr p'^^^ ^n+R'^n were found in the 
final dephased state in the two limiting cases of gi = and g^ = oo. Here we list 
their tails for large R only: 

• When gi = 




R+l 



Cft = { [-^^ ) ' ^^"^ < ^ ' (348) 

, when gf > 1 , 



• When gi = oo 



'i\R 

K2j 



' cos[i2arccos(5/)] , when 5/ < 1 , 

Cr = { ( i\R , , (349) 

I ^2^j , when 3/ > 1 , 



The oscillatory correlation function in the ferromagnetic phase, gf < 1, obtained 
after dephasing from a fully disordered initial state, gi = 00, is a clear qualita- 
tive indication of a non-thermal steady state. Notice similarity of this oscillatory 
correlation function to the correlators (153,171) obtained after a linear quench. 

Dephasing to the GGE instead of a thermal state docs not mean that some 
quantities, which are not very sensitive to fine details of GGE, cannot behave like 
in a thermal state. For instance, as shown in Ref. [179], time correlations of the 
order parameter decay exponentially at the same rate as in a canonical ensemble 
with an effective temperature determined by energy density pumped to the system 
by the sudden quench. 

We can conclude that after a sudden quench in the quantum Ising chain observ- 
ables which are local in fermionic representation relax to GGE. 



3.6. GGE and Luttinger model (LM) 

In this Section, following Ref. [177], we consider local relaxation to GGE in the 
integrable Luttinger model (LM) [180]. This model describes low energy properties 
of a wide class of ID Tomonaga-Luttinger liquids [181]. The Hamiltonian of the 
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LM is 



Hlm = Ho + H2 + Hi, (350) 
Ho = vfpY, ■ ^i(.P)Mp) ■■ , (351) 



H2 = ^Y.92i<l) Mq)JLiq) , (352) 



L 



Hi = \Y.'-'^o.{q)M-q): , (353) 

q,a 

where the Fermi operators satisfy {'i/'a(p)) V'^(p')} — ^p,p'^a,i3 with a, /3 = L,R and 
anticommute otherwise. Anti-periodic boundary conditions ipa{x + L) = —ipa{x), 
where ipa{x) = e'^^^Va (p) / with sr = —sl = 1, result in a non-degenerate 
ground state. The quantised momentum p = 27r(n— l/2)/L with integer n is "half- 
integer". The current operators Jaiq) = Yip '■ i^aip + Q)'^a{p) where q = 27rm/L 
with integer m. 

The currents obey the Kac-Moody algebra [Jaiq), Jfi{(l')] = 2§'^9+9',o'^a/3 which 
allows one to introduce bosonic operators 

boiq) = -ii27T/\q\L)^/^[eiq)jRi-q) - e{-q)JL{q)] , (354) 
hl{q) = i{2T,/\q\L)^l^[e{q)jR{q) - e{-q)JL{-q)] (355) 



for g 7^ 0. Moreover, there are two conserved operators N = Nr + Nl and J = 
Nr — Nl, where = Ja(0). The Hamiltonian i^LM is quadratic in ha, 6q but not 
diagonal. It is diagonalised by the Bogoliubov transformation 

Kq) = bo{q) cosh ip{q) + bl{-q) sinh(g) , (356) 
bHq) = bo{-q) sinh<^(g) + 6J(q) cosh(g) (357) 

with tanh(^(g) = g2{2)/[vF + g4iq)] and becomes 

i^LM = Xl^^(g)kl b\q)b{q) + TT-yjv-^ + ttvj— , (358) 

where v{q) = [{vp + 54(g))^ - dUq)]^^^, vn = v{0)e'^^^^\ vj = t;(0)e-2^(o) _ ^he 
integrals of motion are N , J, and the bosonic occupation numbers b^{q)b{q). 

We consider an interaction quench where both couplings (72(9) and g4{q) are 
suddenly switched on at i = 0. Before the quench Hi,m = Hq is a non-interacting 
Fermi liquid and its ground state is a vacuum for the operators bQ{q). After some 
algebra, see Ref. [177], one obtains a one-body correlation function 

C^^{x,t) = (0| e**^-"Vfl(^)V'i?(0)e-'*^^" |0) . (359) 

It is interesting to compare two asymptotes of the function in the thermodynamic 
limit: 

C^Bix,t) « yo/2t;t)^ ^ when |x| < 2i;(0)t , (360) 
2tt{x + la) 



August 10, 2010 0:9 Advances in Physics revised2 



Advances in Physics 



99 



and 



27r(x + ia) 



when |a;| 3> 2v{0)t 



(361) 



where 7 = sinh 299(0). The correlation function relaxes from the initial Fermi liquid 
form for |x| ^ 2v{0)t to a final non-Fermi liquid form when |x| <C 2v{0)t, but the 
limiting form C^^{x,t 00) has a different exponent 7^ than in the ground state 
of the LM. 

Thus in any bounded region J7 within the quasiparticle horizon, 



the correlation functions are relaxed to a steady state. Given the quadratic form 

of the bosonic Hamiltonian (358) and the general discussion in Section 3.4, it is 
not quite surprising to find that the correlators in the steady state can be obtained 



where X{q) are fixed by initials conditions. We can conclude that any finite subsys- 
tem within the quasiparticle horizon appears relaxed to GGE. 

3.7. GGE and hard-core bosons 

The integrable hard-core bosons in one dimension [182] were realised experimen- 
tally in the seminal experiment [3] , where the momentum distribution in an expand- 
ing cloud of atoms was observed to have a stationary but non-thermal distribution. 
Relaxation in this system was considered in Ref. [174], where it was shown that 
the correct momentum distribution fj. for one-dimensional hard-core bosons can 
be obtained from a GGE for an equivalent fermionic representation of the hard- 
core bosons. However, in a subsequent paper [183] it was demonstrated that the 
stationary state preserves information not only on momentum distribution, but 
also on momentum correlations which are missing in the GGE description. Here 
we outline the argument of Ref. [183]. 

For t < a system of hard-core bosons was in the ground state in a harmonic 
trap potential, 



x\ < 2v{0)t , 



(362) 



from a GGE 




(363) 



Ho = H + Vtrap , 



(364) 



where Vtrap = J dx V{x)p{x) , p{x) is a density operator, and 




(365) 



is the trap potential. At t = the trap potential V{x) is switched off. 




(366) 



and the bosons expand with the translationally invariant Hamiltonian H. 
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It is convenient to make a Jordan- Wigner transformation 



exp 



ITT 



dy p{y) 



(367) 



where the operators '4){x) and ^{x) correspond to spinless fermions and bosons re- 
spectively: {^"(3;), V'Uy)} = Vp{x) 1 iy)] = ^{x — y). The free Hamiltonian becomes 



j_aP_ 

2m dx"^ 



tp{x) , 



but the density operator retains its form 

p{x) = (p^{x)ip{x) = ip^{x)ip{x) . 
Since the fermionic occupation numbers in momentum space 

nk 



i^k'^k , 



(368) 



(369) 



(370) 



where ipk = (27r) J dx e *'^^^(a:), mutually commute and commute with 
H = J dk^rik, they are the integrals of motion to be included in GGE, see the 
discussion in Section 3.3. Neither their expectation values (n^) nor expectation 
values of their products depend on t. In particular, for Snk = rik — (nk) the Wick 
theorem gives a correlator 



{5nk5nk')t 



(371) 



Since the initial ground state in a harmonic trap is not translationally invariant, 
the right hand side is generally non-zero when k' ^ k. 

For a harmonic trap V{x) = with I = (mwo)"^^^ the correlation function 
can be written as 



N-l 



(372) 



n=0 



where (j)n{k) are eigenfunctions of the harmonic oscillator in momentum represen- 
tation. When A/" » 1 we can approximate, see Ref. [183], 



{rik) = ~y^~p"' {Srikdnk') 



siiP [{k - k')R] 
7r2(A; - A;') 



l\2 ' 



(373) 



where kpl = R/l = V2N and ^ kp- There are non-negligible correlators 

between different fermionic occupation numbers. They originate from the broken 
translational invariance of the initial state and are conserved after the opening of 
the trap (366). 

In an experiment like Ref. [3] one can in principle measure bosonic correlation 
functions. Since the Jordan- Wigner transformation between bosons and fermions 
is non-local, the bosonic correlations are in general different than their fermionic 
counterparts. However, Ref. [183] demonstrates that bosonic occupation numbers 
fk tend to the conserved fermionic occupation numbers: fk{t — > 00) = n^. Thus 
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when t ^ oo ell the non-trivial initial correlations between the integrals of motion 
can be measured as bosonic momentum correlations. 

When t oo GGE becomes p*^*-"'^ = AAexp(— J dk Xkfk)- By construction, 
this ensemble "predicts" correct expectation values Trp'^'^^fk = Uk, but it ignores 
any momentum correlations. Thus, in accordance with the discussion in Section 
3.3, GGE cannot be used to predict occupation numbers in momentum space. 
However, the general argument in Section 3.4 proves that, after dephasing for 
t ^ oo, p^^^ = Af exp{— J dk Xkfk) accurately predicts observables within a finite 
fermionic support Q, in real space. 

The local relaxation in Section 3.4 cannot always be taken for granted as demon- 
strated by another example in Rcf. [184]. In that paper, instead of switching off 
the trapping potential, they consider a sudden quench of the trap frequency in Eq. 
(365) from an initial coq to a lower frequency coi: 

(374) 

This problem is solved exactly by a scaling transformation of the initial wave 
function. The wave function is coherently breathing with frequency wi/2. There 
are periodic revivals of the initial state every tt/loi demonstrating the absence of 
dephasing at the revivals. This anomaly originates from the fact that the non- 
interacting Jordan- Wigner fermions in a harmonic potential have commensurate 
eigcncncrgics being integer multiples of ui, so there is no way anything could 
irreversibly dephase when t — )■ oo. Thus one of the basic assumptions in Section 
3.4 is not satisfied and there is no local relaxation to any GGE. However, one 
could argue that a small anharmonicity of the trapping potential would alter this 
conclusion ^. 



t=o 

uo — > . 



3.8. GGE and the Bose-Hubbard model 

The Bosc-Hubbard model (229) is not intcgrabic, but it is close to intcgrability in 
some regimes or limits of parameters. It seems to relax to a non-thermal steady 
state when close to integrability, and thermalise otherwise. 

For instance, Ref. [185] considers a quench from J = to a large J ^ n well in 
the superfluid phase. This quench jumps between two limits of the model where 
it is integrable. The state after the quench relaxes locally to a non-thermal steady 
state. It will be shown below that the steady state is a GGE. This is an example 
that even a non-integrable model can locally relax to a non-thermal steady state 
when it is close to integrability. 

Before the quench, when the Hamiltonian is Hq = ^ J^j ~ the initial 

state is the non-Gaussian Mott ground state |n, n, . . . ) with exactly n bosons at 
each site. Suddenly, at i = a finite hopping rate J is switched on. We assume 
J ^ n so large that the interaction term can be neglected and 

N 

H Ri - J^(a]_^^aj + a]aj+i) = - 2J^ cos(/c)a|.afc (375) 
j=l k 



'-The limiting case of tJi = 0, i.e. tiic opening of tiic trap (366) considered in tiiis Section, results in 
a continuous dispersion k'^ /2rn of the non-interacting .Jordan- Wigner fermions. The integration over the 
continuous k yields irreversible local relaxation when t — > oo. 

This limiting case suggests that even for wi > there may be local relaxation in between the periodic 
revivals such that a local state pii{t) oscillates between the initial pn{0) at the revivals and Tinp'^'^^ 
in between. When oJi — > then the revival period n/wi stretches to infinity allowing enough time for 
irreversible local relaxation. 
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is just the hopping term, where = XljLi o,je~'^^^ j^fN is an annihilation operator 

in pseudomomentum representation, = a^a^ are the integrals of motion of this 
quadratic B.. Since before the quench bosons were localised in the Mott state, their 
conserved occupation numbers are the same for all quasimomenta, Tr p(t) = n. 

We are interested in local obscrvables with support in a finite subsystem of 
the lattice. The Hamiltonian (375) is quadratic as assumed in Section 3.4, but 
the initial Mott state is not Gaussian. Nevertheless, Ref. [185] shows that in the 
thermodynamic limit N ^ oo and t oo (here the order of limits is important) the 
reduced density matrix of a block of 5 sites converges to a product of Gaussian 
states 

5 g-A a]aj 



lim lim pnit) = T] — (376) 



and A = ln(l + 1/n). This limiting steady state can be also obtained by reduction 
of the generalised Gibbs ensemble, 

lim lim pn{t) = Trn±p^^^ , (377) 



where the ensemble is 



^GGE 



Efe A alok p- A a]aj -\ a]aj 



(1 + n)^ (1 + n)^ J-J- 1 + n 

j=i 



This GGE is clearly different from a thermal state oc e"^^*'^'' because uj^ = 
2J(1 — cosfe) is not a constant. We can conclude that, in the integrable limit 
J ^ GO of the non-integrable Bose-Hubbard model, the non-Gaussian initial Mott 
state relaxes locally to a steady non-thermal GGE. 

A similar problem was addressed in Ref. [186] by numerical simulations of small 
lattices up to 12 sites. Figure 37 shows probability distributions Pa in the diagonal 
ensemble p in Eq. (313). A quench from the Mott-insulator to the superfiuid phase 
prepares a clearly non-thermal diagonal ensemble, see Fig. 37d, but a small quench 
within the superfiuid phase prepares a diagonal ensemble that looks more like a 
canonical ensemble, see Figs. 37a,c. 

A sudden quench in the opposite direction, i.e. from the superfiuid to Mott- 
insulator phase, was studied by exact diagonalisation [187], the adaptive time- 
dependent density matrix renormalisation group [169], and Bogoliubov theory 
[121]. The model is trivially integrable at J = 0, where unitary evolution fac- 
torises into a product over lattice sites W- e^*"^ ^"^^^)/^. Since the quartic interac- 
tion nj{nj — l)/2 is an integer for any Fock state, the evolution is periodic in time 
with a period of 2tt. Any initial wave function revives periodically and there is no 
question of any relaxation for t — > oo. However, a non-zero hopping rate J > 
leads to relaxation of these periodic revivals or even ovcrdamped relaxation, see 
Ref. [128], Section 2.19.3 and the 2D experimental results in Fig. 29. The numer- 
ical data collected in Ref. [187] suggest that there are two distinct regimes: when 
J is close to the Mott-supcrfluid transition, then for simple observables the final 
steady state is a thermal state with an effective temperature determined by energy 
pumped into the system by the sudden quench, but when J is perturbatively small, 
then the steady state is not thermal and it retains memory of the initial state. The 
last result is further supported by Fig. 37b from Ref. [186]. 
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Figure 37. Here Ui and Uf are initial and final interaction strengths in the Bose-Hubbard Hamiltonian 
(229) in a sudden quench Ui Uf. The critical point between the Mott and superfluid phases is Uc = 3.3. 
The plots show probabilities pn in a mixed time-averaged density matrix after a quench, p = pn\n){n\, 
versus the excitation energy LOn — Eq of an eigenstate \n). Panels a and c show roughly exponential 
dependence of p„ on energy, like in a canonical ensemble, while panels b and d indicate strong memory of 
the initial conditions. The results come from numerical calculations on a periodic lattice of 8, . . . , 12 sites 
with a density of 1 particle per site. (Figure from Ref. [186]) 



The qualitatively different behaviour in the two regimes may be explained by a 
simple model in Ref. [121]. In this model, excitations of the Mott insulator state 
|1, 1, 1, . . . ) are particles, created by p\ and located at doubly occupied sites, and 
holes, created by h\ and located at empty sites. When expanded to second order 
in the operators p, h the Bose-Hubbard Hamiltonian (229) becomes a quadratic 
Hamiltonian which can be diagonalised by a Bogoliubov transformation to a sum 
of non-interacting Bogoliubov quasiparticles, H2 = Ylik a^kl\ dlk.cx- This quadratic 
model, which neglects any interactions between quasiparticles 7^,0-, is accurate for 
small J. Thus the general argument in Section 3.4 predicts local relaxation to a 
non-thermal GGE with the integrals of motion 71^^^ = l\ alk,a- 

In this model, thermalisation can occur due to quasiparticle interactions. For 
instance, quartic terms of the form '^\^k+q/2l-k+q/2lo could make quasiparticle 
population equilibrate. Deep in the Mott regime, where the Mott gap ~ 1 in the 
quasiparticle spectrum uoj. is large as compared to the quasiparticle bandwidth 
~ J, this process is not compatible with energy conservation, but it becomes in- 
creasingly effective as J approaches the gapless critical point and it may explain 
the thermalisation observed close to criticality. However, a limited accuracy of the 
leading quadratic model was discussed already in Ref. [121] and it is also not known 
how these qualitative observations depend on a finite system size. 

3.9. GGE and a system solvable by Bethe Ansatz 

Reference [188] considers a one-dimensional fermionic lattice Hamiltonian ^ 

H= -^(^4^^cs + li.c}j +V^nsns+i (379) 

s s 



-'^ Transport was analysed in the system (379) in the first paper in Ref. [189]. 
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Figure 38. A time-averaged momentum distribution for the final V = 10, i.e., momentum distribution in 
the time-averaged diagonal ensemble p, and thermal momentum distribution in a canonical ensemble. The 
inset shows a time dependence of (n^) together with a (dotted) thermal average. (Figure from Ref. [188]) 



with nearest neighbour interaction V at half filling. This Hamiltonian is not 
quadratic and cannot be mapped to any quadratic Hamiltonian, but it is exactly 
solvable by Bethe ansatz. The model has a quantum critical point at 14 = 2 sepa- 
rating a Luttinger liquid (when V < 2) from a charge density wave insulator (when 
V > 2). Ref. [188] considers open chains up to 100 sites pushed out of equilibrium 
by sudden quenches from Vq to V and evolved in time using the Lanchos time 
evolution method [190] and the adaptive DMRG [169]. The quantity of interest is 
momentum distribution after time t: 

1 ^ 

i^k) = J^Y. ^''^"""^ Tr p{t) clcn . (380) 

m,n=l 

A general conclusion drawn from the numerical simulations in Ref. [188] is that 
the momentum distribution relaxes to a steady state, but the steady state is not a 
canonical ensemble, see Fig. 38. 

Since the Hamiltonian (379) cannot be mapped to any quadratic Hamiltonian, 
Section 3.4 cannot tell us what might be the small set of integrals of motion la in 
the GGE (322). In principle, any quantum system has as many integrals of motion 
as the dimension of its Hilbert space given by the projectors la = \oi){oi\ on the 
eigenstates \a) of its Hamiltonian. However, GGE build out of these projectors is 
equal to the time-averaged diagonal ensemble 

/O = Y^Pa \(y){a\ (381) 

a 

which may be not a tractable statistical description of the state. 

Another universal set of integrals of motion are integer powers of a finite Hamilto- 
nian: Im = with m = 1, . . . , dim(ff). This set is as huge as the set of eigenstate 
projectors but, at least in the quench considered in Ref. (379), it can be accurately 
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truncated to a small number of M leading powers: 

pGGE ^ ^exp I - > ^ Xrr.H'^ 1 (382) 




with A's fixed by the conserved moments Trp'^^^if™' = Trp(0)-fr'". The ensemble 
(382) can be rewritten as p'^^^ = where 



p{E) = AAexp (- X^eA 

\ m=l J 



(383) 



This form reveals that the essence of the generalised GGE (382) simply is the 
assumption that pa in Eq. (381) is a smooth localized function of the energy E^, 
which can be accurately approximated by e.g. an exponent of a short polynomial in 
Ea- Good quality of this approximation for reasonably small M, with M = 1 being 
a canonical ensemble, is demonstrated in Fig. 39 which shows p{E) multiplied by 
density of eigenstates p{E), i.e., P{E) =p{E)p{E). 



3.10. Eigenstate thermalisation hypothesis (ETH) 

Since a state after a sudden quench often has a narrow energy distribution in a 
final Hamiltonian, see e.g. Fig. 39, and the energy distribution p{E) is conserved, 
it may be more accurate to describe a thermal state of an isolated quantum system 
by a microcanonical ensemble 

P(^) = A^li'f l.^-^^^K^^' (384) 
I U, otherwise , 

Here {H) is average energy after a quench and /S.E is a finite energy window. In this 
framework thermalisation of an observable O means that the exact time-averaged 
diagonal ensemble p = 'Yl,aPaW){(A ^-^^ the microcanonical ensemble (384) give 
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the same expectation value of O: 



''^^ Pa{oe\0\a) = M (a|0|a) , 



(385) 




where \a) are eigenstates of the Hamiltonian. 

One can imagine many different scenarios that might in principle imply the 
equality (385), but there is only one that does not depend on the initial state, 
except that p^'s have to be localised around the average energy {H). The eigenstate 
thermalisation hypothesis suggested in Ref. [163] proposes that the expectation 
value (a|0|a:) of a few-body observable O in an eigenstate |q;) with energy Ea of a 
large many-body Hamiltonian is the same for all eigenstates in the energy window 



ETH clearly implies the equality (385) independently of the details of Pa, pro- 
vided that it is only localised around a given average energy. This is what is ex- 
pected from thermalisation: the relaxed state should not depend on the initial state, 
except for its average energy. 

The thermalisation mechanism implied by ETH is very simple: any eigenstate in 
the energy window gives the same expectation value of a few body observable O 
or, even stronger, any single eigenstate in the window provides the same average as 
the microcanonical ensemble. This mechanism is very different from the ergodicity 
required of a classical system. In a non-intcgrablc isolated quantum system, each 
eigenstate of the Hamiltonian implicitly contains a thermal state. The coherence 
between the eigenstates initially hides it, but then time evolution reveals it after 
dephasing. 

At present, there are no general theoretical arguments supporting ETH, but 
there are some results for restricted classes of Hamiltonians: Deutsch in Ref. [163] 
showed that ETH holds for an integrable Hamiltonian perturbed by a matrix from 
a random Gaussian ensemble, nuclear shell calculations have shown that individual 
wave functions reproduce thermodynamic predictions [191], some quantum systems 
whose classical versions are chaotic satisfy ETH in the semiclassical limit [192]. 
More generally, ETH follows from the Berry conjecture [86, 163] believed to hold 
in such systems [193]. 

In order to see if thermalisation occurs in a generic isolated quantum system 
and, in particular, if it occurs due to ETH, Ref. [194] considers a non-integrable 
system of 5 hard core bosons on a finite lattice of 21 sites in Fig. 40a described by 
the Hamiltonian 



with the nearest neighbour repulsion strength U = O.IJ. The bosons are initially 
confined to the bottom-right part of the lattice in the ground state of a confined 
Hamiltonian. At t = they are allowed to tunnel to the initially empty top-left 
part. This set-up is a quantum analogue of an inflated balloon pierced inside a 
vacuum chamber to see that the released air will soon uniformly fill the chamber 
and velocity distribution of air molecules will relax to the Maxwell velocity distri- 
bution. Reference [194] considers momentum distribution. For instance. Fig. 40b 
shows relaxation of the distribution at fc,,. = to the microcanonical ensemble. 
Fig. 40c shows that the momentum distribution in the dephased diagonal ensem- 
ble is clearly different from the initial momentum distribution, indistinguishable 
from the microcanonical momentum distribution, but significantly different from 



Ea ± AE. 




(386) 
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Figure 40. Panel A shows the lattice in the Hamiltonian (386). The bosons are initially prepared in the 
ground state of the sub-lattice in the bottom-right corner and then, at t = 0, released through the marked 
link. Panel B shows how the momentum distribution n{kx = 0) relaxes to the value in the diagonal ensemble 
which is the same as predicted by the microcanonical ensemble. Here n{kx, ky) = 'Y^^ ^ £—^f'(ri—rj)/L 
with lattice size L = 5 and n{kx) = X]fe n{kx,ky). Panel C, compares the initial momentum distribu- 
tion n{kx) with the steady state diagonal ensemble which is indistinguishable from the microcanonical 
distribution but significantly different from the canonical one. (Figure from Ref. [194]) 



the canonical one. 

Figure 41 compares results for the non-integrable system of 5 hard-core bosons 
on the lattice in Fig. 40a, which was considered so far, and a similar integrable 
system of 5 hard-core bosons in a one-dimensional chain of 21 sites. The bosons 
were initially prepared in the ground state of the 8-sites at one of the chain's 
ends, and then released by opening the link connecting the end to the rest of the 
chain. Results collected in Fig. 41 show that while ETH holds in the non-integrable 
system, it does not hold in the integrable one. 

Similar results were also obtained in Ref. [195] for interaction quenches in a finite 
one-dimensional chain with hard-core bosons whose integrability can be gradually 
broken by turning on the next-nearest-neighbour hopping and repulsion. When it 
is not integrable the system thermalises and ETH holds, but as the next-nearest- 
neighbour terms are turned off and the system tends to the integrable limit both 
the thermalisation and ETH break down. The transition between the two regimes 
seems to be a smooth crossover, but it is not clear how much this conclusion 
depends on the finite system size. 
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Figure 41. Left panels (A,B,C) correspond to the non-integrable system, and the right panels {D,E,F) 
to the integrable one. Panels A and D compare momentum distribution n{kx) with the microcanonical 
distribution, and with distributions in two typical eigenstates close to the average energy. The four distri- 
butions collapse in the non-integrable case (ETH holds), but they differ significantly in the integrable case 
(ETH does not hold). The top parts of panels B and E show expectation values of n{kx = 0) in eigenstates 
at different energies. The non-integrable case (B) looks like a smooth curve, but the integrable case (E) 
shows large fluctuations which do not resemble a smooth curve. The lower parts of panels B and E show 
that energy distributions are very similar in both non-integrable and integrable case. Finally, the upper 
parts of panels C and F show a focus on expectation values of n{kx = 0) in the eigenstates |a) close to 
the mean energy, and the lower parts of panels C and F show probabilities pa in these eigenstates. (Figure 
from Ref. [194]) 



3.11. Dynamics of relaxation to a steady state 

So far we have discussed mainly the nature of the asymptotic steady state. Finally, 
it is time for the dynamics of relaxation to this state. In this Section we briefly 
review two References [196, 197] where the relaxation process was studied in some 
detail. 

In the first of them [196] they study a one-dimensional fermionic Hubbard model 
H = Yl ^^^-4^,. + UJ2 Ut - I) Ui - I) ■ (387) 

ija i 

Here the hopping matrix Vij corresponds to a semi-elliptic density of states 
p{e) = \/4 — e^/ (27r). The system is initially in the ground state of the non- 
interacting Hamiltonian, U = for i < 0, and then at t = the Coulomb repulsion 
is switched to a finite U. The relaxation following this sudden quench was studied 
in Ref. [196] within the time-dependent dynamical mean-field theory. Evolution of 
two observables, the double occupancy d{t) = Trp{t)ni^nii and the Fermi surface 
discontinuity An(t), is shown in Fig. 42. For both U ^ 1 and U <^ 1, the two 
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Figure 42. Double occupancy d(t) and Fermi surface discontinuity An(t) after quenches to [7 < 3 (left 
panels) and U > 3.5 (right panels). The horizontal dashed lines in panel c indicate the quasi-stationary 
value predicted in Ref. [143], and the horizontal arrows in panel a mark corresponding thermal values. 
(Figure from Ref. [196]) 



quantities pass through a plateau where, as predicted in Ref. [143], the system pre- 
thermahses in a quasi-steady state, before it finahy relaxes to a thermal state. The 
size of the plateau depends on U. More detailed analysis shows that it is minimal 
at ?7 ~ 3.2 marking a dynamical phase transition between the weak and strong 
coupling regimes where the thermalisation is postponed until after the prether- 
malization on intermediate time-scales. As remarked in Ref. [196], it is not clear 
whether and how this phenomenon is related to the existence of an equilibrium 
thermodynamic phase transition at Uc = 4.76, but the relaxation is the fastest at 
the crossover between the strong and weak coupling regimes. 

A similar problem was considered in Ref. [197] in the antiferromagnetic XXZ 
model (solvable by Bethe ansatz) and the antiferromagnetic XZ model (integrable 
by mapping to a quadratic Hamiltonian) , 



s 

Hxz = (2<^.\i + Aa,V,\i) , 



(388) 



(389) 



They studied unitary time evolution of the antiferromagnetic order starting from 
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Figure 43. Relaxation time and oscillation period as functions of anisotropy A in the XXZ and XZ models. 
In both models the relaxation time is the shortest at the critical point. (Figure from Ref. [197]) 

the highly non-equihbrium Neel state. The order vanishes exponentially with relax- 
ation being oscillatory or non-oscillatory, depending on the anisotropy parameter 
A. As demonstrated by the data collected in Fig. 43, the relaxation is the fastest 
near the quantum critical point contrary to the usual notion of critical slowing 
down. This effect can be explained by the gapped excitation spectrum away from 
the critical point. The relaxation is dominated by scattering between high-energy 
excitations introduced to the system through the highly excited initial state. The 
gap restricts the phase space available for scattering making relaxation slower. This 
leads to increasing relaxation time with increasing gap and minimal relaxation time 
near the critical point where the gap is zero. 

In conclusion, in all models reviewed in this Section relaxation is the fastest either 
at a critical point or close to it. This is where either vanishing or small energy gap 
does not inhibit relaxation. 



3.12. Summary 

In this part we focused on the nature of the apparent stationary states long time 
after a sudden quench with a special emphasis on quadratic models. Dynamics of 
relaxation to the steady state was only touched in Section 3.11. This is a non-trivial 
problem, especially for a large sudden quench, where not only low energy states 
are excited, but also high energy ones which are not described by any universal low 
energy effective theory. This problem should in general be less serious after a linear 
quench where, despite some similarity to a sudden quench in the adiabatic-impulse- 
adiabatic approximation, only low energy modes are excited up to a cut-off set by 
the KZ length ^, and this length itself can be accurately obtained from a universal 
low energy theory. However, sudden quenches are interesting in their own right. 
After the pioneering theoretical work in the 1970's [198], they became recently the 
subject of dedicated experiments on oscillations and dephasing of the superfluid 
state after a sudden quench deep into the Mott insulator phase [2, 199]. 
There are integrable models where the time evolution can be solved exactly, like 



the XY chain [83, 178, 198], ID hard-core bosons [174, 183], or the ^ Hubbard 
chain [175], all of which are diagonalizable by the Jordan- Wigner transformation. 
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but this type of integrability may result in very specific relaxation dynamics [200]. 
In models integrable by the Bethe ansatz it was not possible to extract dynamics in 
general, but with the notable exceptions of the Richardson and Lieb-Linigcr models 
[201]. In view of these difficulties, it is necessary to develop approximate and/or 
numerical methods such as solutions of the dynamics of field theoretical models at 
the renormalisation group fixed points [167, 177, 202], semiclassical theories [203], 
exact diagonalisation [186, 194, 204], timc-dcpcndcnt dynamical renormalisation 
group (tDMRG) in one dimension [78, 169, 187, 188, 204, 205], dynamical mean 
field theory in infinite dimensions [142, 196, 206], and other approximate methods 
[128]. R(K-cntly, new numerical schemes were proposed [207] which go beyond the 
standard tDMRG by focusing from the beginning on the dynamics of the interesting 
local observable rather than a quantum state as a whole. Some of the more universal 
results in the above papers concerning the final relaxed steady states were reviewed 
in this part, but the less universal dynamics of relaxation has been only touched. 

As explained in the introductory Section 3.1, it makes sense to define a mixed 
steady state poo even for an isolated quantum system in a pure state, provided that 
the statistical description is used only for observables O that are sufficiently local 
or coarse-grained. This observation encourages one to apply the usual concepts of 
thermodynamics to pure quantum states. The discussion in Section 3.1 revealed 
the central role played by the diagonal part p = '^aPa\'^)i'^\ °f density matrix 
in the eigenbasis |q) of the time- independent Hamiltonian H after the quench. 
Reference [208] goes further and argues that also for a time-dependent Hamiltonian 
H{t) it is only the diagonal part of the density matrix p(t) = \xjj(t)){ip{t)\ in the 
instantaneous eigenbasis of H{t) that is relevant for thermodynamics, because any 
realistic measurement of a thermodynamic observable O actually measures its time 
average O = TrpO. The time average dephascs to zero any contribution from the 
off-diagonal elements of the density matrix p{t), and for sufficiently coarse-grained 
observables the dephasing may be fast enough to justify using the time-dependent 
diagonal ensemble as a statistical description of the actual underlying pure 
state \tp{t)). This observation motivates definition of the diagonal entropy 



which is in general non-zero, in contrast to the vanishing von Neumann entropy 
S = —TTp{t) log2 p{t) of the underlying pure state. The diagonal entropy does not 
change in adiabatic processes when transitions between instantaneous energy levels 
are negligible and the diagonal probabilities pa do not change. In Rcf. [208] it is 
shown that when the system is initially in thermal equilibrium at temperature T, 
then Sd can only increase or stay the same, in accordance with the second law 
of thermodynamics. This and other properties [208] make the diagonal entropy a 
plausible candidate for the entropy of an isolated quantum system. 

Another Ref. [209] derives microscopic expression for the heat generated in an 
arbitrary process as the energy generated due to transitions between different in- 
stantaneous energy levels. When the initial density matrix is diagonal in the Hamil- 
tonian eigenbasis, then the expression for heat contains only squares of the absolute 
values of the evolution operator which can be interpreted as transition probabilities 
between instantaneous energy levels. It can be shown [209] that when the initial 
density matrix is passive, i.e. its diagonal elements Pa decrease with increasing 
cigcncnergy, then the heat is non-negative. 

A closely related issue is considered in Ref. [210], where a probability distribution 
P{W) is derived for a work W done in an instantaneous quench starting at the 




(390) 



a 
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critical point of a quantum critical system. For instance, in the quantum Ising 
chain, a small sudden quench from the critical transverse magnetic field ^ = 1 to 
a nearby g = 1 + Sg does a work W with the Gamma distribution 



This distribution exhibits an interesting edge singularity with the exponent 
{5g/2'KY depending on the size 5g of the quench. Certainly, a lot of interesting 
research could be done along these lines and those of Ref. [209] especially in the 
context of the Jarzynski theorem [211]. 

Another Reference [212] considers spin susceptibility in a sudden quench of the 
transverse field g in the quantum Ising chain. After the field is instantaneously 
switched from an initial gi to a final gf the transverse magnetisation relaxes locally 
to a stationary value m{gi,gf) = Trp(7|. One can define susceptibility to the initial 
field as Xi = dm/dgi. Exact calculations in Ref. [212] show, among other things, 
that close to the criticality, when g'j fa 1, the quench susceptibility has a logarithmic 
divergence 



where F{gf) is independent of gi. Notice that the divergent susceptibility probes 
not only the ground state, which is known to undergo drastic changes at a quantum 
critical point, but also the excited states. 

These are but a few examples which did not quite fit into the main story of 
Section 3, but which demonstrate that the subject is far from closed and we can 
expect exciting new developments in the near future. 

4. Open problems 

The main open problems appear to be: 

• Once the non-adiabaticity of quantum critical systems has been well established, 
we must at last face the problem how to prevent the unwanted excitation in the 
adiabatic quantum state preparation. The non-linear quench or the inhomoge- 
ncous transition in Sections 2.5 and 2.8 respectively are the first proposals in 
this direction, but we certainly still need more ideas and more work to overcome 
this problem. 

• Due to the quasiparticle horizon effect, the entropy of entanglement in the excited 

state grows linearly with time and makes simulations with classical computers 
a formidable task. On the other hand, the excited high energy eigenstates ren- 
der low energy effective theories inaccurate. Thus an accurate solution of the 
relaxation dynamics in a large non-integrable system remains a challenge. 

These and other open problems may require another review in the near future. 
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